1. Introduction

It is hoped that x-ray Phase Contrast Imaging (XPCi) will provide a generational improvement in the effectiveness of mammography [1]. To our knowledge, the only in vivo mammography program is in progress in Trieste, Italy, using the SYRMEP beam line [2]. This program has provided mammograms of improved spatial resolution and detail visibility compared with conventional mammography. It cannot, however, be considered a viable tool for clinical screening due to its reliance on a synchrotron source.

An alternative XPCi technique employing laboratory sources was suggested by Olivo et. al [3, 4] in 2007. This technique is known as coded aperture XPCi and has since been under continuous development within the radiation physics group at UCL (see references [5, 6, 3] for example). This technique has been demonstrated experimentally and validated theoretically in the aforementioned references. We are now building a pre-prototype coded aperture XPCi system in order to demonstrate the efficacy of the technique using in vitro human breast tissue samples. In order to design the system and verify the experiments, it is necessary to model the entire imaging system, including the interaction of x-rays with tissue. The small refractive index contrast of tissue combined with the unpolarised x-ray source mean that a full electromagnetic calculation for the scattered x-rays can be avoided. Furthermore, the short wavelength of x-rays relative to typical cell structure dimensions means that a geometrical optics model is often sufficient. This is important as a rigorous scalar calculation of the scattered field would require prohibitively large computational resources. In this paper, we thus attempt to establish conditions under which a geometrical optics approximation can be employed to model a coded aperture XPCi system.

As with models of other, related imaging systems [7, 8, 9, 10, 11], we consider phase objects whose optical thickness is at least piece–wise smooth. Work was done by Peterzol et. al [12] to determine the limits of validity of the geometrical optics approximation for free space propagation type XPCi systems. Such an analysis has not been performed for coded aperture XPCi systems. The link between geometrical and wave optics is by no means a new area of research. Keller [13] was the first to show that geometrical optics need not be limited to modelling objects with smoothly varying refractive index. He showed that geometrical optics is an approximation to wave optics which can be made more accurate by the inclusion of higher order terms. A good account of this technique is given by James [14]. In this paper we calculate higher order terms to show how the geometrical optics and wave optics solutions vary in predicting coded aperture XPCi images. The paper is arranged as follows. We first present the wave optics model and show how it can be implemented efficiently. We then derive the geometric optics model before showing how a source of finite width can be introduced into the system. We then apply the developed theory to the particular example of an infinite cylinder. By employing the stationary phase approximation to the diffraction integrals which result from the wave optics model, we derive the geometrical optics model, thus showing how the two models are related. Finally we show some numerical examples and show conditions under which the geometrical optics model may be accurately employed.

2. Wave optics model

We consider first the wave optics model of the imaging system depicted in Fig. 1. Normally a sample would be placed on the detector side of the sample apertures however we initially consider the sample free case. Following the method employed by Olivo and Speller [16] we use Fresnel-Kirchhoff diffraction theory to calculate the field incident upon the detector apertures. We consider initially a single point source at position (x_s,0,−z_so) emitting a spherical wave at wavelength λ. Previous experiments have shown [16] that modelling the system at the source’s average energy gives a good prediction of the image. The assumption of a point source will be relaxed later. Assuming the exp(−iωt) sign convention, the field at position P = (x,y, z_od) may be given by [16]:

where

and 𝓐 represents the transmitting regions of the sample apertures. In addition, x, ξ and z are defined in Fig. 1 and (ξ, ψ, z) and (x,y, z) form right handed coordinate systems. The integration over ψ can be performed by noting that the apertures have no dependence upon ψ. We must thus evaluate:

evaluating this integral with limits at ±∞ contradicts the Fresnel approximation used to obtain Eq. (1) and should be solved using the theory of distributions [15]. This problem can however be avoided by noting that the kernel of the integral in Eq. (3) is a rapidly oscillating function which lends itself to asymptotic evaluation by the method of stationary phase. According to the method of stationary phase [14, Pgs. 29–34], an integral of the form:

 

Fig. 1. Schematic diagram of imaging system including reference frames used in the paper. Note that (x̄,ȳ,z), (ξ, ψ, z) and (x, y, z) all form right handed coordinate systems. The imaging system is assumed to have no y dependence. Note that dL is defined by the displacement between the detector apertures and the projection of the sample apertures onto the detector apertures.

Download Full Size | PDF

where g(x) has a single first order stationary point, x ₀, such that g′(x ₀) = 0, g″(x ₀) ≠ 0, can be approximated as:

in the limit of large k. Applying this approximation to Eq. (3) we find that

the role of this term is to ensure energy conservation and give the incident field the correct phase relationship with y. This result is also obtainable using Fourier theory applied to distributions [15] which reveals that Eq. (6) is in fact the solution to Eq. (3) [17]. This enables us to write Eq. (1) as:

where we now introduce the periodic function T(ξ) to represent the transmission function of the sample aperture. It is now easy to include the effect of a phase object with phase function ϕ(ξ) by following an approach similar to that of Arfelli et. al [18]. The total field at the detector apertures may be found according to:

where ℬ is the extent of the object.

3. Efficient evaluation of wave optics field

We now turn our attention to how the expression in Eq. (7) may be efficiently evaluated. As T(ξ) is a periodic function with period L, it can be represented as a complex Fourier series written in general as:

which upon substitution into Eq. (7) yields

As also suggested by Engelhardt et. al [10], the Fast Fourier Transform (FFT) can be used to efficiently evaluate this expression. In particular, starting with the definition of the discrete Fourier transform [19]:

by allowing x′z_so/(z_so + z_od) to take on values κL/(2N) the summation in Eq. (10) may be evaluated, for a finite number of terms, by constructing a vector of the form:

where

and finally taking the Fourier transform of the vector in Eq. (12). Noting also that the coefficients C_n may also be evaluated using the FFT, Eq. (10) may be evaluated very efficiently.

The second term in Eq. (8) must, in general, be evaluated numerically unless the object has a phase function permitting analytic evaluation. It was found that Gaussian Quadrature integration [19] provided accurate results.

4. Geometrical optics model

Olivo and Speller [6, 20] have previously used geometrical optics to model the coded aperture XPCi system. Their approach used a “forward” technique where photons emitted by the source were traced through the system. Photons could be blocked by an aperture, refracted by a sample or both. The number of photons reaching a particular pixel represent the signal detected by that pixel. We now consider the ray optics approach in a more formal manner in order to relate it to the wave optics approach. For the remainder of this section we consider only non-trivial rays which are transmitted by the sample aperture. We consider here a first order geometrical optics. It has been shown by Keller [13] and later by James [14] that geometrical optics may be extended to include higher order terms which represent what is usually termed diffraction. Here we consider only the first order terms of the geometrical optics approximation. The trajectory of a light ray is described by the expression [21]:

where r is the position vector of a point on the ray, s the length of the ray, n the refractive index of the medium and 𝓢 defines a wave front of constant phase, ie, 𝓢 (ξ) = constant. It is evident from this that we assume rays are deflected in the ξ direction only. Consider a phase object as depicted in Fig. 2. We define the phase function, ϕ(ξ), as

where n(ξ, z) is the refractive index at position (ξ, z) and we have assumed that rays make only small angles, θ_i, with the z-axis. The angle by which the ray is deflected in then given by:

With reference to Fig. 1, we can say that a ray emitted at angle θ_i to the z-axis will intercept the ξ-axis at position ξ = z_so tan(θ_i) and, if deflected by an object, will intercept the x-axis at position

The phase of the ray at the detector apertures is calculated by taking into account the phase introduced by the object and the distance travelled in free space according to the Fresnel approximation. The amplitude of the ray must be such that energy is conserved. In particular, the time average power propagating in a small pencil of rays emanating from the source must remain constant. The ratio between ray amplitudes at z = z_od and z = 0 is thus given by:

 

Fig. 2. Diagram illustrating the phase function ϕ(ξ) of a phase object of extent z_ob in the z direction.

Download Full Size | PDF

5. Modelling a finite size source

Secs. (2)–(4) show how to calculate the field incident upon the detector apertures. The detected signal is found by integrating the intensity of x-rays transmitted by the detector apertures and incident upon a particular pixel. In general, the pth transmitting region of the detector apertures is given by [p LM−LM/4+dL,pLM+LM/4+dL] where dL is the displacement of the detector apertures relative to the projection of the sample apertures as shown in Fig. 1 and M = (z_od + z_so)/z_so is the system magnification. We assume that the pixels are aligned as shown in Fig. 1 such that a single pixel entirely covers a single transmitting region of the detector apertures. Before calculating the signal detected by each pixel, we introduce a source of finite size in the x̄ direction. The brightness is described by P(x̄) which we will take to have a Gaussian profile. We can then take the signal of the pth pixel to be given by:

where, for mathematical convenience we have assumed that the source brightness profile limits the effective source size rather than the limits of integration. By making the substitution P(x̄) = exp(-(x̄/σ)²), Eq. (19) may be expressed as:

where

and erf(z) is the error function

Equation 20 shows that K(x) may effectively be considered as a pixel sensitivity function. Figure 3 shows plots of K(x) for a variety of source Full Widths at Half Maximum (FWHM). This shows how a broad source leads to a broad K(x) thus diminishing the sensitivity of the system to fine variations in the intensity caused by phase variations in the object.

 

Fig. 3. Plot of K(x) for some values of source FWHM, given in the legend in μm. Values of M = 1.25 and L = 40μm were used.

Download Full Size | PDF

6. XPCi model of a dielectric fibre

6.1. Wave optics model

We now apply the results of previous sections to model the XPCi image of a cylindrical fibre. This problem has been considered previously by Olivo and Speller [6] in order to verify experimental results. We consider a non-absorbing cylinder of radius R, refractive index n = 1−δ, parallel to the ψ-axis centered upon (ξ, z) = (ξ ₀,0). Absorbing materials can be modelled by writing n = 1 − δ + iβ thus introducing an attenuation term in Eq. (8). We have opted to set β to 0 to simplify the following analysis. Note that δ is of the order of 10⁻⁶ to 10⁻⁷ for the range of x-ray energies and materials which we consider here. The phase function, ϕ(ξ), may thus be calculated as:

We consider first the wave optics calculation. We must evaluate the second term of Eq. (8) after substituting Eq. (23) into it. The bounds of integration are found by taking the intersection of the transmitting part of the sample aperture and the cylinder. Without loss of generality we consider three cases depicted in Fig. 4 where the transmitting part of the sample aperture is assumed to be centered upon ξ = 0. Note that the cases depicted in Fig. 4 do not limit the cylinder radius, all that is important is where the cylinder boundaries lie relative to the transmitting regions of the apertures. A cylinder covering more than one sample aperture could be modelled using a combination of the cases depicted in Fig. 4. In practice our system employs a series of apertures to simultaneously image a wide field of view. For clarity, we consider here a single sample/detector aperture pair and scan the object to obtain its image. Images obtained in this way will be equivalent to those obtained in practice only when photons are not scattered between differing pre-sample/detector aperture pairs. Only a simple extension is required to model the practical system as is shown at the end of Sec. (6.2). An analysis of when this approximation is valid is given in Sec. (6.3).

 

Fig. 4. Three cases which must be considered in order to evaluate the integral in Eq. (8).

Download Full Size | PDF

The integration may be evaluated numerically however it is instructive to analytically evaluate by approximation. We start by writing Eq. (8) as

where

where Ω = [ξ ₀ − R, ξ ₀ +R] ∩ [−Lη/2,Lη)/2] and η is the fill factor of the sample apertures. We now attempt to find asymptotic solutions, for large k, to the integrals in Eq. (25) by again applying the stationary phase approximation. In this case we must also consider the end points of the integrals. It is shown by James [14, Pgs. 29–34] that when the integration in Eq. (4) is evaluated over the interval [a, b] and g′(x) is non-zero and finite at the end points, a term

must be included. It is shown by Murray [22, Pg. 77] that if one of the end points is a stationary point, then Eq. (26) should be omitted and a factor of 1/2 included in the term of Eq. (5). We now turn our attention to evaluating leading terms in the asymptotic expansions of the integrals in Eq. (25). We start by defining the functions g ₁ (ξ) and g ₂ (ξ) and finding their derivatives as:

it is easy to verify that when M, δ and z_od are limited to those values experienced in practice, g′₁ (ξ _1,0) = 0 has a unique solution for every value of x′. This solution must in general be calculated numerically. This may be done efficiently by evaluating γ = [γ_i] = [g′₁(ξ_i)+x′/z_od] where ξ_i = [ξ_i] is a discretisation of the domain [ξ ₀ − R + ε, ξ ₀ + R − ε] for some small ε. This corresponds to case (2) in Fig. 4 where the entire cylinder is illuminated and thus rays are refracted to all values of x′. In cases (1) and (3), the bounds of integration are affected by the sample apertures which in turn affects the values of x′ to which rays are refracted. The stationary point of g ₁ for a particular x′ may be found by interpolation with γ as the abscissa. This enables the leading term in the expansion to be calculated as

examination of g′₂ shows that g ₂ has a single stationary point at ξ _2,0 = x′ /M. In the case that x′/M is within the bounds of integration of U ₂, the following term is contributed by the stationary point ξ _2,0:

supposing then that Ω = [a, b] and that neither a or b are stationary points of g ₂, the next term in the asymptotic expansion may be found as

which, in the special case where b = −a, becomes

6.2. Relationship between wave and geometrical optics models

We show here how the wave and geometrical optics solutions are related. As θ_i and α in Eq. (17) are small we can write Eq. (17) as

where Eq. (23) has been substituted into Eq. (16) to find α. This expression is identical to g′₁ = 0 in Eq. (27), the stationary phase condition for integral U ₁. It is then easy to verify that assuming identical incident field conditions, substituting Eqs. (23) and (16) into Eq. (18) results in the same magnitude as CΓ(y)I _1,0. Furthermore, substitution of the phase contributions from the phase object and the Fresnel approximation for free space propagation result in the same phase as in CΓ(y)I _1,0. This shows that the leading term in the asymptotic expansion of U ₁ gives the same field as the geometrical optics approximation to the refracted field. Examination of CΓ(y)l _2,0 = 1/(z_so + Z_od) exp (ik((y ² + (x − x_s)²)/(2(z_so + z_od)) + z_so + z_od)) shows that this is the geometrical optics field of the light which reaches the detector apertures without being refracted by the cylinder or blocked by the sample apertures. Closer examination of CΓ(y)I _2,1 shows that this is the field due to diffraction at the edges of Ω. Note that this quantity becomes infinite at the edges of the geometrical projection of Ω onto the detector apertures. This non-physical result can be remedied by modifying the stationary phase solution [23] however this is beyond the scope of this work.

Table 6.2 shows how the intensity and complex amplitude, for the geometrical and wave optics models respectively, are calculated in each region defined in Fig. 5. Note that the geometrical optics solution provides the intensity of the field whilst the wave optics solution provides the complex amplitude of the field.

 

Fig. 5. Diagram illustrating the three regions which must be considered when analysing the field incident upon the detector apertures.

Download Full Size | PDF

Table 1. Summary of components contributing to the field and intensity in the different regions for the case of wave and geometrical optics solutions respectively. Note that objects illuminated by adjacent sample apertures may be modelled by replacing each term I _1,0 and I _1,0 ² with summations, Σ_i I ⁱ _1,0 and Σ_i(I ⁱ _1,0)², over all objects i respectively.

View Table

6.3. Examples and analysis

The validity of the presented model and technique depend on the tendency for photons to be scattered between adjacent sample/detector aperture pairs. It is however possible to develop a minimum bound upon the separation of apertures required to maintain validity. If a cylinder of radius R is placed with its centre at ξ = 0 in the imaging system of Fig. 1, its edge will be projected onto the position x′ = MR in the space of the detector. We are interested in knowing how quickly the field scattered by the cylinder decays away from x′ = MR. Assuming that the edge of the cylinder is illuminated, photons are refracted to values of x′ approaching ∞ and are described by the term I _1,0 defined in Eq. (28). Photons reaching a position x′ ≫ MR must be incident upon the cylinder for a value of ξ very close to, but not exceeding R. By writing ξ = R − ε, ε > 0, in Eqs. (27) it is easy to find a simple analytic expression giving I _1,0 for x ≫ MR as ε tends to 0. It is then simple to show that I _1,0 ² will reduce by two orders of magnitude at a position x′ = RM + Δx′ where:

Figure 6 shows contours of Δx′ for values of R and δ encountered in practice. Δx′ may be considered the minimum separation of adjacent sample/detector aperture pairs to ensure detector apertures principally detect photons originating from their associated sample aperture. The above analysis considers only a point source. A source of finite width may be considered by noting the definition of x′ in Eqs. (2) and thus adding (W/2)z_od/z_so to Δx′, where W is the detector FWHM.

 

Fig. 6. Contours of Δx′ for a variety of values of R and δ. This diagram effectively shows the minimum separation required between adjacent sample/detector aperture pairs to ensure detector apertures principally detect photons originating from their associated sample aperture.

Download Full Size | PDF

Previous studies [6] have shown that coded aperture XPCi contrast is increased by reducing the fraction of detector pixel exposed to directly incident radiation. This however leads to an increase in the exposure time as fewer photons reach the pixel. In this work we have thus chosen a displacement, dL, equal to half of the transmitting width of the detector apertures, thus exposing half of the pixel to directly incident radiation. We used a sample aperture periodicity of L = 40μm along with z_so = 1.6m and z_od = .4m to match the dimensions of an experimental system currently under construction. The simulations were performed for a photon energy of 100keV.

Figure 8 shows the intensity incident upon the detector apertures as calculated by the wave optics and geometrical optics solutions for a point source illuminating a cylinder. The cylinder has a value of δ = 10⁻⁷, a radius of 5μm and was situated with its axis at ξ = − 5μm. As is expected, the wave optics intensity exhibits oscillations resulting from interference between different field components. The geometrical optics solution is physically impossible as the sharp edge occurring at x = 0 would require the field to contain infinite spatial frequencies. Consideration of the angular spectrum of a propagating aperiodic field shows that such a field would require evanescent waves which, in our case, would have negligible magnitude such a distance from the sample apertures.

Figure 9 compares the directly calculated wave optics intensity to that calculated using the stationary phase approximation. As explained in Sec. (6.2), singularity anomalies arise in this solution which have been neglected. This plot shows that apart from these anomalies, the approximate solution agrees well with the directly calculated intensity. One can use the components which comprise the approximate solution to determine when the geometrical optics and wave optics solutions converge. This is however made difficult by the singularity anomalies present in the approximate solution and so we have opted to use a more pragmatic approach as outlined below.

One can envisage that when a source of finite width is employed, the XPCi signals predicted by the geometrical and wave optics models should converge. There are two explanations for why this should be the case. The first explanation observes that the point source intensity incident upon the detector apertures is convolved with the magnified source profile which in our case is Gaussian. This is equivalent to applying a low pass filter to the intensity distribution causing the oscillations in the wave optics intensity and the sharp transition in the geometrical optics intensity to be smoothed. The second explanation considers that, as shown in Sec. (5), a source of finite width may be modelled by a system employing a point source and equivalent detector apertures which cause the pixels to have a spatially dependent sensitivity as described by Eq. (21). Figure 3 shows that as the source broadens, so does the width of the equivalent detector aperture sensitivity function. Because of energy conservation, one would expect the geometrical and wave optics XPCi signals to converge as the source broadens. In particular, consider the plot shown in Fig. 7. This shows the difference between the intensities, incident upon the detector apertures, predicted by wave and geometrical optics. This signal has a zero mean value as required by conservation of energy. The coded aperture XPCi signal thus depends upon the domain over which the field intensity is integrated by the detector pixel. As the sensitive part of each detector pixel increases, or equivalently, as the source broadens, the geometrical and wave optics signals thus tend to converge.

 

Fig. 7. Plot of the difference between intensities calculated using the wave optics (full expression evaluated numerically) and geometrical optics approximations. The sensitive region of the pixel is shaded. Simulation parameters were the same as in Fig. 8.

Download Full Size | PDF

Previous studies [6] have shown that the maximum XPCi signal for a cylindrical object using the system described in this section occurs when the cylinder is positioned at approximately ξ ₀ = −R. This is demonstrated in Fig. 10 where wave and geometrical optics signal traces have been plotted for a cylinder of radius 5μm and δ = 10⁻⁶. The signals have been normalised by the signal for the object free case. These plots demonstrate how the signal traces converge as the FWHM of the source increases. It also shows how the peak of each trace is in the vicinity of ξ ₀ = −R, as expected. Simulations run over a range of radii, values of δ and source FWHM show that the peak of the signal trace does indeed occur in the region of ξ ₀ = −R. This is suggests a good way of assessing the difference between the wave and geometrical optics XPCi signals as the two signals are likely to vary most at the peak. We thus calculate an error term, ε(−R), where ε(ξ ₀) = |I _WO(ξ ₀)/I _WO ^N − I _GO ^N(ξ ₀)/I _GO ^N|, and I _WO(ξ ₀) and I _GO(ξ ₀) are the XPCi signals for the geometrical and wave optics (full expression evaluated numerically) cases respectively, for a cylinder at position ξ ₀. I _GO ^N and I _WO ^N are the object free XPCi signals for the geometrical and wave optics cases respectively.

 

Fig. 8. Plot of the intensity of the field incident upon the detector apertures for the geometrical and wave optics (full expression evaluated numerically) solutions. Simulation parameters used were R = 5μm and ε = 10⁻⁷, all other parameters were as described in Sec. (6.3).

Download Full Size | PDF

 

Fig. 9. Plot of the intensity of the field incident upon the detector apertures as calculated using the exact and approximate wave optics formulations. Simulation parameters were the same as in Fig. 8.

Download Full Size | PDF

Before proceeding to calculate ε it is useful to note that some approximations can provide further insight into the problem. In the case of ξ ₀ = −R, g ₁ in Eq. (27) can be well approximated by $g_1\approx -\mathrm{x\prime }\xi /z_{\mathrm{od}}-2\delta \sqrt{-2\mathrm{\xi R}}$ for x′ > −MR, but not too close to −MR. This approximate form leads to a solution of ξ _1,0 = −2δ ² Rz_od ²/x′² for the stationary point of g ₁. Substitution of ξ _1,0 back into the approximate forms of g ₁ and g″₁ show that both of these functions have a dependence upon δ ² R rather than each of these independently. This suggests that it is reasonable to expect ε for a particular source FWHM to be constant for constant values of δ ² R. This is indeed the case as was verified by a large number of simulations, a small selection of which are shown in Fig. 11. This significantly simplifies the task of determining the source size for which the geometrical and wave optics signals converge. Figure 12 is a contour plot of e as a function of source FWHM and δ ² R. The important conclusion which we can draw from this is that for our particular choice of z_od and z_so, as we expect a source to have a FWHM of around 50μm, the geometrical optics model will provide results consistent with those of the wave optics model. This result will make it feasible to model much larger objects.

 

Fig. 10. Plot wave optics (WO, full expression evaluated numerically) and geometrical optics (GO) XPCi signal traces for a cylinder of radius 5μm and δ = 10⁻⁶, for three different values of source FWHM. The cylinder is scanned from ξ ₀ = −L/4 − R to ξ ₀ = 0. The signals have been normalised by the signal for the object free case.

Download Full Size | PDF

 

Fig. 11. Plots of ε against δ ² R for three different values of δ and a point source (left) and a source of FWHM 50μm.

Download Full Size | PDF

 

Fig. 12. Contour plot of the error between the normalised XPCi signals as calculated by geometrical and wave optics models. Source FWHM is on the vertical axis and δ ² R is on the horizontal axis.

Download Full Size | PDF

7. Conclusions

In this paper we have outlined the two most widely used techniques for modelling XPCi systems: wave and geometrical optics. We have used the theory developed to model the image of an infinite cylinder in a coded aperture XPCi system. This problem has practical significance as it can be tested experimentally. For this particular problem, we show how the geometric and wave optics models are related. we then show how this theory can be used to develop a guide for when the two techniques can be trusted to give consistent results.