Phase retrieval using coherent imaging systems with linear transfer functions

Optics Communications 234 (2004) 87–105 www.elsevier.com/locate/optcom

Phase retrieval using coherent imaging systems with linear transfer functions David Paganin a,*, Timur E. Gureyev b, Konstantin M. Pavlov a, Robert A. Lewis a, Marcus Kitchen a a

b

School of Physics and Materials Engineering, Monash University, Vic. 3800, Australia Commonwealth Scientific and Industrial Research Organisation, Manufacturing and Infrastructure Technology, Private Bag 33, Clayton South, Vic. 3169, Australia Received 27 June 2003; received in revised form 23 January 2004; accepted 6 February 2004 D.P. would like to dedicate this paper to the memory of Armando Paganin

Abstract We consider the problem of quantitative phase retrieval from images obtained using a coherent shift-invariant linear imaging system whose associated transfer function (i.e., the Fourier transform of the complex point-spread function) is well approximated by a linear function of spatial frequency. This linear approximation to the transfer function is applicable when the spread of spatial frequencies, in a two-dimensional complex wavefield, is sufficiently narrow when compared to the characteristic length of variation of the transfer function for an imaging system taking such a wavefield as input. We give several algorithms for reconstructing both the phase and amplitude of a given two-dimensional coherent wavefield, given as input data one or more images of such a wavefield which may be formed by different states of the imaging system. When an object to be imaged consists of a single material, or of a single material embedded in a substrate of constant thickness, the phase-amplitude reconstruction can be performed using a single image. As a first application of these ideas, we write down an algorithm for using a single diffraction-enhanced image (DEI) to obtain a quantitative reconstruction of the projected thickness of a single-material sample which is embedded within a substrate of approximately constant thickness. This algorithm is used to quantitatively map inclusions in a breast phantom, from a single synchrotron DEI image of the same. In particular, the reconstructed images quantitatively represent the projected thickness in the bulk of the sample, in contrast to raw DEI images which greatly emphasise sharp edges (high spatial frequencies). Lastly, we point out that the methods presented here are also applicable to the quantitative analysis of differential interference contrast (DIC) images, obtained using both visible-light and X-ray microscopy. Ó 2004 Elsevier B.V. All rights reserved. PACS: 03.65.Vf; 03.65.Wj; 41.60.Ap; 42.30.Rx; 42.30.Wb Keywords: Phase contrast imaging; Phase retrieval; X-ray imaging

*

Corresponding author. Tel.: +61-3-9905-3633; fax: +61-3-9905-3637. E-mail address: [email protected] (D. Paganin).

0030-4018/$ - see front matter Ó 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.optcom.2004.02.015

88

D. Paganin et al. / Optics Communications 234 (2004) 87–105

1. Introduction At optical and higher frequencies, the extreme rapidity of wavefield oscillations implies that one is only able to directly measure the square modulus of the complex amplitude of such a wavefield, averaged over many cycles [1]. To combat this limitation, many means have been developed for so-called ‘‘phase-contrast imaging’’, which may be defined as any means of imaging which is sensitive to phase variations across a given electromagnetic or matter wavefield. These methods for phase-contrast imaging are necessarily indirect, for they measure the modulus of a suitably modified wavefield in order to make qualitative and/or quantitative inferences about the phase of the unmodified wavefield. The most famous method for phase contrast imaging, namely interferometry, is one of many methods for rendering visible the transverse phase shifts in a coherent wavefield, such as might be imprinted upon a plane wave traversing a sample. Other methods include Zernike phase contrast [2], Lorentz microscopy [3], in-line holography [4], differential interference contrast (DIC) microscopy [5], Schlieren imaging [6], focus variation methods [7] and diffraction enhanced imaging (DEI) [8]. After briefly examining the features generic to all such means of phase-contrast imaging (Section 2), we discuss the case of images, obtained using a coherent shift-invariant linear imaging system, for which the transfer function may be taken to be a linear function of spatial frequency (Section 3). This linear approximation to the transfer function (which we term a ‘‘linear transfer function’’) is applicable when the spread of spatial frequencies, in a given two-dimensional complex wavefield, is sufficiently narrow (with respect to the characteristic length of variation of the transfer function of an imaging system taking such a wavefield as input). Note that the assumption of a linear transfer function, while rather strong at first sight, is applicable to a number of realistic imaging systems, as shall be outlined later in the paper. We develop several algorithms for phase-amplitude reconstruction of a given wavefield; these algorithms require, as input, one or more images of the disturbance which are taken using a coherent imaging system possessing a linear transfer function (Section 4). The four-image reconstruction algorithm, which requires as input data four different images of the same two-dimensional wavefield which are obtained using four different states of the imaging system, makes no further assumptions regarding the imaging system beyond those listed previously. However, we show that one can perform the phase retrieval using fewer images, if one makes additional assumptions. In particular, one can make do with a single image provided that: (i) one is using plane waves to image an object which is composed of a single-material sample of interest, embedded in a substrate of approximately uniform thickness; (ii) the projection approximation is valid, i.e., the phase and amplitude shifts of the probe radiation, upon passage through the object, are the same as would be obtained if the objectÕs potential distribution were projected onto a plane perpendicular to the optic axis. Within its restricted domain of validity, the application of such a single-image algorithm is simpler than methods which require two or more images as data. As a detailed case study and experimental implementation of aspects of the general formalism developed in Section 4, Section 5 considers the problem of quantitative diffraction-enhanced imaging. This technique, which has generated much recent interest in the X-ray imaging community, achieves phase contrast by diffracting a given wavefield from a so-called analyzer crystal inclined at an angle to an incident wave-front [8] (see Fig. 1). The physical mechanism for this means of phase contrast, the origins of which date back at least as far as 1959 [9], is the extreme sensitivity of the crystalÕs reflection coefficient to phase gradients in the incident wave [10]. We use this technique to demonstrate quantitative phase-amplitude reconstruction, using a single DEI image of inclusions in a breast phantom. Significantly, by suitable analysis of the phasecontrast effects in a given DEI image, the method is able to generate the absorption radiographs which would have been obtained if (i) the radiation dose was higher and (ii) all scattered radiation was eliminated. The stability of this algorithm is studied, leading to the conclusion that it is robust in the presence of realistic levels of noise in the data. Since it requires only a single image as input, the single-image algorithm

D. Paganin et al. / Optics Communications 234 (2004) 87–105

89

Fig. 1. Schematic of experimental setup for diffraction-enhanced imaging. Monochromatic planar X-ray radiation passes through a sample which consists of a substance of interest embedded in a matrix (which may be vacuum) of approximately uniform thickness. Such a sample is termed ‘‘quasi-homogeneous’’ in the present paper; if the matrix is vacuum, then the object is said to be ‘‘homogeneous’’. Reflection of the transmitted radiation from an appropriately oriented analyzer crystal rejects the wide-angle scattered radiation, and produces image contrast in the coherently scattered component, using the mechanism discussed in the text. After reflection from the analyzer crystal, the DEI image is recorded using a position-sensitive detector.

may be applied to real-time image data obtained using DEI imaging (Section 5) and DIC microscopy (Section 6).

2. General remarks on phase-contrast imaging The most well-known means of coherent phase-contrast imaging is interferometry, which renders visible the transverse phase variations in a given two-dimensional wavefield by first modifying the wavefield through superposing it with a reference beam, and then measuring the resulting intensity. In the language of Fourier optics, and for the case of a plane reference wave, the action of an interferometer may be viewed as greatly amplifying a single spatial frequency in the two-dimensional Fourier transform of the monochromatic wavefield of interest (i.e., multiplying this Fourier transform by 1 þ adðkx
kx0 ; ky
ky0 Þ, where a is a complex parameter, ðkx ; ky Þ are the Fourier space coordinates dual to the real space coordinates ðx; yÞ over the plane of interest, kx0 ; ky0 are the transverse spatial frequencies of the reference plane wave expðikx0 x þ iky0 yÞ over the plane of interest, and d is the Dirac delta). Other means of phase-contrast imaging include in-line holography, Zernike phase contrast, Schlieren imaging and DIC. What all of these methods have in common is this: like interferometry, they suitably modify the two-dimensional Fourier spectrum of a wavefield of interest and then measure the square modulus of the modified field, in order to render visible the previously invisible effects of the transverse phase variations of the unmodified wavefield. Regarding the four examples just cited: (a) in-line holography renders phase visible via free space propagation through a distance z, which is equivalent to multiplying the two-dimensional Fourier transform of theffi unpropagated complex disturbance by the phase factor (transfer qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi function) expðikz zÞ ¼ expðiz k 2
kx2
ky2 Þ, where k  2p=k and k is the radiation wavelength [11]; (b) Zernike phase contrast renders phase visible by multiplying the lowest spatial frequency components of the two-dimensional Fourier transform of a disturbance by a complex constant of modulus unity, this being the effect of the central phase-shifting ‘‘spot’’ in the simplest incarnation of the Zernike phase plate [2]; (c) Schlieren imaging renders phase visible by using a ‘‘knife edge’’ to set to zero all spatial frequencies, in the two-dimensional Fourier transform of a wavefield disturbance, over a given half-plane whose edge passes through the origin of coordinates in Fourier space; (d) Differential interference contrast renders phase visible by interfering a given wave with a phase-biased and slightly displaced copy of itself [5] (see Section 6): using the Fourier shift theorem, one can easily show that this is equivalent to multiplying the Fourier transform of a disturbance by 1 þ expðið/0
kx dx
ky dyÞÞ, where /0 is the phase bias and ðdx; dyÞ are Cartesian coordinates of the displacement vector of the original wavefield from the phase-biased copy of itself. One may classify the various methods of phase-contrast imaging as either ‘‘qualitative’’ or ‘‘quantitative’’. The former class, utilised by most practitioners of phase contrast imaging, is primarily interested in

90

D. Paganin et al. / Optics Communications 234 (2004) 87–105

qualitative (e.g., morphological) characteristics of the phase-contrast image produced by a given sample in a given imaging system. The latter class of imaging, commonly termed ‘‘phase retrieval’’ or ‘‘phase reconstruction’’, seeks to infer quantitative phase and amplitude information from a given set of phasecontrast images; if desired, this information may subsequently be related to the quantitative structural properties of a sample with which the phase-retrieved wavefield has interacted. 3. Coherent image formation using a shift-invariant linear imaging with a linear transfer function Let us consider the problem of imaging a coherent complex scalar wavefield, the spatial part of which is denoted by wðx; y; zÞ, where x; y; z are Cartesian coordinates in three-dimensional space. We suppress the harmonic time dependence expð
ixtÞ of the coherent wavefield, where x denotes angular frequency and t is time. The plane z ¼ z0 is chosen to be the entrance surface of a given imaging system; this may either be a physical surface or a mathematical one, and is defined to be the plane over which the wavefield to be imaged assumes its ‘‘input’’ value. Denote by wIN ðx; yÞ  wðx; y; z0 Þ a given complex disturbance over the entrance surface of our coherent imaging system. The output two-dimensional wavefield wOUT ðx; yÞ, produced by the imaging system A, may be written in operational terms as ^ A wIN ðx; yÞ; wOUT ðx; yÞ ¼ O

ð1Þ

^ A is an operator which completely characterises the imaging system. Such an imaging system is said where O ð1Þ ð2Þ to be linear [12] if a linear combination of two inputs, say wIN ðx; yÞ and wIN ðx; yÞ, leads to the same linear combination of corresponding outputs h i ^ A awð1Þ ðx; yÞ þ bwð2Þ ðx; yÞ ¼ aO ^ A wð1Þ ðx; yÞ þ bO ^ A wð2Þ ðx; yÞ; O ð2Þ IN IN IN IN where a and b are arbitrary complex numbers. If, furthermore, a given coherent linear imaging system is shift invariant [12], then a transverse spatial shift in input yields a shift in the corresponding output, with the shift in the output being proportional to the shift in the input. The action of an arbitrary coherent shift-invariant linear imaging system may be represented by a convolution integral [12], allowing us to rewrite Eq. (1) as wOUT ðx; yÞ ¼ wIN ðx; yÞ  Kðx; yÞ;

ð3Þ

where  denotes two dimensional convolution, and Kðx; yÞ is the convolution kernel (complex point-spread function) which characterizes the coherent shift-invariant linear imaging system. By the convolution theorem of Fourier analysis [12], one may Fourier transform Eq. (3) with respect to x and y, to yield the product ~ ~ ~ w OUT ðkx ; ky Þ ¼ wIN ðkx ; ky ÞKðkx ; ky Þ:

ð4Þ

Here, a tilde over R R a function w denotes its two-dimensional Fourier transform using the convention ~ x ; ky Þ dkx dky , and ðkx ; ky Þ are the previously introduced Fourier-space wðx; yÞ ¼ ð2pÞ
2 expðikx x þ iky yÞwðk ~ x ; ky Þ is termed the ‘‘transfer function’’, or the ‘‘propagator’’, of coordinates dual to ðx; yÞ. The quantity Kðk the imaging system. ~ ðkx ; ky Þ, of the input two-dimensional coherent Assume that the two-dimensional Fourier transform w IN wavefield, is non-negligible over a patch of Fourier space which is sufficiently small that one may safely ~ x ; ky Þ in Eq. (4), this linearization being with respect to the Fourier-space coordinates kx ; ky . In linearize Kðk regions of Fourier space for which this linearization yields an inaccurate approximation to the transfer ~ ðkx ; ky Þ be sufficiently small for this error in the linearization to be irrelevant: it function, we require that w IN ~ ðkx ; ky Þ is small, since it is the ~ x ; ky Þ is poor in regions where w does not matter if the estimate for Kðk IN ~ ~ product wIN ðkx ; ky ÞKðkx ; ky Þ which appears in (4). Thus (4) becomes

D. Paganin et al. / Optics Communications 234 (2004) 87–105

h i 0 0 0 0 ~ ~ ~ w ðk ; k Þ  w ðk ; k Þ Kðk ; k Þ þ ðk
k ÞN þ ðk
k ÞX ; x y OUT x y IN x y x y x y

91

ð5Þ

where we have assumed that the previously mentioned patch of reciprocal space is centred about ðkx0 ; ky0 Þ, ~ x ; ky Þ with its first-order Taylor series approximation about this point. In (5), we allowing us to replace Kðk have introduced the symbols     o ~ o ~ Kðkx ; ky Þ Kðkx ; ky Þ N ; X ð6Þ ok ok 0 0 0 0 x

y

ðkx ;ky Þ¼ðkx ;ky Þ

ðkx ;ky Þ¼ðkx ;ky Þ

for the kx and ky coordinates of the two-dimensional Fourier-space gradient rk of the transfer function ~ x ; ky Þ at ðk 0 ; k 0 Þ (see Fig. 2). We emphasize that this process, of linearizing the transfer function by Kðk x y writing it as a linear function of spatial frequency, is a central (and rather restrictive) approximation of the present paper. Note, for future reference, that: (i) the transfer function in the square brackets of (5) will be termed a linear transfer function; (ii) a coherent shift-invariant linear imaging system, for which (5) is a good approximation, will be said to have a linear transfer function. ~ can be made to As is evident from the geometric construction of Fig. 2, one of the coordinates of rk K vanish, at the point ðkx0 ; ky0 Þ, by a simultaneous rotation of both the real-space coordinates ðx; yÞ and the Fourier space coordinates ðkx ; ky Þ, about their respective origins, by an anti-clockwise angle h ¼ tan
1 ðX=NÞ. Let the rotated real-space coordinates be denoted by ðX ; Yffi Þ, with the rotated reciprocal pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi space coordinates given by ðn; gÞ. Further, N becomes T  signðNÞ N2 þ X2 in the new coordinate system, while X vanishes. Thus (5) becomes h i ~ ~ ðn; gÞ KðA ~ w ðn; gÞ  w ; A Þ þ ðn
A ÞT ; ð7Þ n g n OUT IN where the point ðkx0 ; ky0 Þ in the old Fourier-space coordinate system is ðAn ; Ag Þ in the rotated system. Take the inverse Fourier transform of Eq. (7) with respect to n and g, to give wOUT ðX ; Y Þ  SwIN ðX ; Y Þ
iT

o w ðX ; Y Þ; oX IN

~ n ; Ag Þ
An T ; S  KðA

ð8Þ

where use has been made of the Fourier derivative theorem [12], using the following convention for the twodimensional Fourier transform

ky

η

Ω

k'y Αη O

Ξ

∼ ∇kΚ

θ ξ

Αξ θ

k'x

kx

Fig. 2. The Fourier space dual to the plane ðx; yÞ, which is transverse to the optic axis z, is coordinatized by the ordered pair ðkx ; ky Þ. ~ x ; ky Þ in Eq. (4), is linearized about The transfer function (propagator) for a shift-invariant coherent linear imaging system, namely Kðk  ~ denote the gradient of Kðk ~ x ; ky Þ at ðkx0 ; ky0 Þ i.e., rk K ~  ^kx oKðk ~ x ; ky Þ=okx  ~ x ; ky Þ= the point ðkx0 ; ky0 Þ. Let rk K þ ^ky oKðk ðkx ;ky Þ¼ðkx0 ;ky0 Þ oky jðkx ;ky Þ¼ðkx0 ;ky0 Þ , where ^kx and ^ky respectively denote unit vectors in the kx and ky directions. In the ðkx ; ky Þ coordinate system, ~ ¼ ðN; XÞ makes an angle h with the positive kx axis. Rotate the coordinate system by an anti-clockwise angle h ¼ tan
1 ðX=NÞ to rk K ~ ¼ ðT ; 0Þ has a zero projection along the g axis, where give a newpffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi coordinate system ðn; gÞ, which is such that rk K ffi T ¼ signðNÞ N2 þ X2 . In passing from the old to the new coordinate system, ðkx0 ; ky0 Þ ! ðAn ; Ag Þ. We can then write the linearized transfer function as it appears in the square brackets of Eq. (7).

92

D. Paganin et al. / Optics Communications 234 (2004) 87–105

wðX ; Y Þ ¼ ð2pÞ


2

Z Z

~ gÞ dn dg: expðiX n þ iY gÞwðn;

ð9Þ

Write bothpffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi wIN ðX ; Yffi Þ and wOUT ðX ; Y Þ in terms of their intensity and phase, i.e., let 2 wIN ðX ; Y Þ ¼ IIN ðX ; Y Þ expðiuIN ðX ; Y ÞÞ, where IIN ðX ; Y Þ  jwIN ðX ; Y Þj and uIN ðX ; Y Þ  arg wIN ðX ; Y Þ; IOUT ðX ; Y Þ and uOUT ðX ; Y Þ are similarly defined. Eq. (8) then becomes   pffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffi i o o loge IIN þ T uIN ; IOUT exp ½iðuOUT
uIN Þ ¼ IIN S
T ð10Þ 2 oX oX where the ðX ; Y Þ dependence has been suppressed for clarity. Take the square modulus of (10), to yield ( 2 ) ðST  Þ oXo uIN
ImðST  Þ oXo iloge IIN jS j þ 2Re h  2  2 IOUT ¼ IIN ; þjT j2 oXo uIN þ 14 oXo loge IIN

ð11Þ

where Re ðST  Þ and Im ðST  Þ respectively denote the real and imaginary parts of ST  , and a superscript asterisk denotes complex conjugation. On account of the derivatives with respect to X , the image IOUT will evidently show relatively strong ‘‘derivative’’ contrast for phase and intensity variations along the X direction, with relatively less contrast for identical fluctuations along the Y direction. This completes our account of the forward problem, of calculating intensity images output by a given coherent shift-invariant linear imaging system, which possesses a transfer function which may be approximated as a linear function of spatial frequency. We now turn to the associated inverse problem, of deducing the phase and amplitude/intensity of the coherent disturbance input into the said imaging system, given one or more intensity images which are output by the system.

4. Quantitative phase-amplitude reconstruction Here, we show how one may use Eq. (11) to recover a given input complex wavefield wIN ðx; yÞ, given as data one or more images of this wavefield which are obtained using a coherent linear imaging system whose action is well approximated by Eq. (5). For brevity, this system will henceforth be called the ‘‘imaging system’’. 4.1. Reconstruction using four images Suppose that our imaging system is used to form four different images of the same input wavefield pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi wIN ðX ; Y Þ ¼ IIN ðX ; Y Þ expðiuIN ðX ; Y ÞÞ. Subject to the approximations used in deriving Eq. (11) , a given state of the imaging system is completely characterised by the pair of complex numbers ðS; T Þ. Let the four ð1Þ ð2Þ ð3Þ ð4Þ images IOUT ; IOUT ; IOUT ; IOUT be formed using four different states of the imaging system, which are characterised by the known complex ordered pairs ðSj ; Tj Þ, j ¼ 1; 2; 3; 4. Note that each of these states must have the same X axis, a condition which applies to all multiple-image reconstruction algorithms presented in this paper. Use (11) to write down the four equations:    2  2 o o ðjÞ       u
ImðSj Tj Þ loge IIN þ Tj C ; j ¼ 1; 2; 3; 4; IOUT ¼ IIN Sj þ 2Re ðSj Tj Þ ð12Þ oX IN oX h i 2 2 where C  ðouIN =oX Þ þ 14 ðo loge IIN =oX Þ . Exclude C from the j ¼ 1; 2 cases of (12); also, exclude C from the from the j ¼ 2; 3 cases and from the j ¼ 3; 4 cases. This leaves the three equations:

D. Paganin et al. / Optics Communications 234 (2004) 87–105

93

b

2Re aij oXo uIN
Im aij oXo loge IIN þ IINij þ cij ¼ 0; ði; jÞ ¼ ð1; 2Þ; ð2; 3Þ; ð3; 4Þ; 2 ðjÞ ðiÞ I I jSj j S jSi j2 c 
: aij  TSii
Tjj ; bij  OUT2
jOUT 2 ; 2 2 ij Ti j jTi j jTjj jTj j

ð13Þ

Eqs. (13) may be considered as three simultaneous linear equations in the three unknowns ouIN =oX , o loge IIN =oX and 1=IIN , which can be solved for these three quantities to give: o 1 aI23 b34 c12 þ aI12 b23 c34 þ aI34 b12 c23
aI12 b34 c23
aI34 b23 c12
aI23 b12 c34 uIN ¼ ; I R I R I R I R I R I oX 2 aR 12 a34 b23 þ a23 a12 b34 þ a34 a23 b12
a12 a23 b34
a23 a34 b12
a34 a12 b23

ð14aÞ

I R I R I R I R I R I aR 23 a34 b12 þ a34 a12 b23 þ a12 a23 b34
a34 a23 b12
a12 a34 b23
a23 a12 b34 ; I R I R I R I R I R I aR 12 a34 c23 þ a23 a12 c34 þ a34 a23 c12
a12 a23 c34
a23 a34 c12
a34 a12 c23

ð14bÞ

IIN ¼

R R R R o aR b c þ aR 34 b12 c23 þ a23 b34 c12
a12 b34 c23
a23 b12 c34
a34 b23 c12 loge IIN ¼ R12 I23 34 ; I R I R I R I R I oX a12 a34 b23 þ aR 23 a12 b34 þ a34 a23 b12
a12 a23 b34
a23 a34 b12
a34 a12 b23

ð14cÞ

I where aR ij  Re aij and aij  Im aij . The two non-degeneracy conditions for Eqs. (14a) and (14c) are Tj 6¼ 0  and b  ða  a Þ 6¼ 0, where b  ðb12 ; b23 ; b34 Þ and a  ða12 ; a23 ; a34 Þ. Similarly, the non-degeneracy conditions for Eq. (14b) are Tj 6¼ 0 and c  ða  a Þ 6¼ 0, where c  ðc12 ; c23 ; c34 Þ. Eq. (14a) , which yields an image influenced only by the phase of the input wavefield, can be integrated with respect to X to yield the reconstructed input phase uIN ðX ; Y Þ. Eq. (14b), and the exponential of the X integral of (14c), can be used to obtain independent estimates for the reconstructed input intensity IIN ðX ; Y Þ. Having reconstructed both IINffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðX ; Y Þ ffi and uIN ðX ; Y Þ, one can reconstruct the complex disturbance wIN ðX ; Y Þ ¼ p IIN ðX ; Y Þ expðiuIN ðX ; Y ÞÞ which was incident upon the imaging system. We close this section by pointing out that an alternative means of reconstructing IIN ðX ; Y Þ is to multiply the left side of (14b) by (14c), and then add the result to (14b). The resulting expression is ð1 þ o=oX ÞIIN ¼ GðX ; Y Þ, the right hand side of which is a known function G of ð1Þ ð2Þ ð3Þ ð4Þ S1 ; T1 ; S2 ; T2 ; S3 ; T3 ; S4 ; T4 ; IOUT ; IOUT ; IOUT and IOUT . This equation may be stably solved using the Fourier
1 derivative theorem to give IIN ¼ F
1 ½ð1 þ inÞ FG, where F denotes Fourier transformation with respect to X, and F
1 denotes the corresponding inverse transformation (cf. Section 4.4.2).

4.2. Reconstruction using three images Suppose that one uses the linear-transfer-function imaging system in three different states ðSj ; Tj Þ, ð1Þ ð2Þ ð3Þ 2 j ¼ 1; 2; 3, to form the three images IOUT ; IOUT ; IOUT . Further, suppose that the jTj j C term in (12) may be neglected, which corresponds to either the input complex amplitude or the transfer function being slowly varying functions. One then has the three equations:    2 o o ðjÞ     IOUT ¼ IIN Sj þ 2Re ðSj Tj Þ u
ImðSj Tj Þ loge IIN ; j ¼ 1; 2; 3: ð15Þ oX IN oX Exclude IIN from the j ¼ 1; 2 pair of equations, and then exclude IIN from the j ¼ 2; 3 pair of equations, to give 2Re -ij oXo uIN
Im -ij oXo loge IIN þ fij ¼ 0; ði; jÞ ¼ ð1; 2Þ; ð2; 3Þ;  2 ðiÞ ðjÞ ðiÞ ðjÞ -ij  IOUT Sj Tj
IOUT Si Ti ; fij  IOUT Sj 
IOUT jSi j2 :

ð16Þ

These two equations are linear in the unknowns ouIN =oX and o loge IIN =oX , and may be solved for these quantities to yield:

94

D. Paganin et al. / Optics Communications 234 (2004) 87–105

o f12 Im -23
f23 Im -12 u ¼ ; oX IN 2½ðIm -12 ÞðRe -23 Þ
ðRe -12 ÞðIm -23 Þ

ð17aÞ

o f12 Re -23
f23 Re -12 loge IIN ¼ : oX ðIm -12 ÞðRe -23 Þ
ðRe -12 ÞðIm -23 Þ

ð17bÞ

The non-degeneracy condition is x  x 6¼ 0, where x  ðx12 ; x23 Þ. Using Eq. (17), one may easily obtain the intensity IIN ðX ; Y Þ, phase uIN ðX ; Y Þ and complex disturbance wIN ðX ; Y Þ of the wavefield incident on the imaging system; these expressions will not be written down here. 4.3. Reconstruction using two images Suppose that one uses the imaging system in two different states ðSj ; Tj Þ, j ¼ 1; 2, to form the two images ð1Þ ð2Þ IOUT and IOUT . Further, suppose that both the jTj j2 C term and the ImðSj Tj Þo loge IIN =oX term in (12) can be neglected; such an approximation is made in existing quantitative analyses of diffraction enhanced images [8], a point we shall return to in Section 5.1. The above approximation amounts to assuming that (i) the 2 input complex amplitude or the transfer function are sufficiently slowly varying that jTj j C may be neglected in (12); and (ii) o loge IIN =oX is not so slowly varying that IIN may be considered a constant (if IIN is constant, then we have the ‘‘pure phase’’ case dealt with in Section 4.4.1). With these approximations, one writes down the two equations    2 o ðjÞ    u ; j ¼ 1; 2; IOUT ¼ IIN Sj þ 2Re ðSj Tj Þ ð18Þ oX IN which are readily solved for IIN ðX ; Y Þ and ouIN ðX ; Y Þ=oX : ð1Þ

ð2Þ

o 1 I jS2 j2
IOUT jS1 j2 uIN ¼  ð2Þ OUT ; ð1Þ oX 2 IOUT Re ðS1 T1 Þ
IOUT Re ðS2 T2 Þ ð2Þ

IIN ¼

ð19aÞ

ð1Þ

IOUT Re ðS1 T1 Þ
IOUT Re ðS2 T2 Þ jS2 j2 Re ðS1 T1 Þ
jS1 j2 Re ðS2 T2 Þ

:

ð19bÞ

Here, the non-degeneracy condition is ð1Þ

IOUT ð2Þ

IOUT

6¼

Re ðS1 T1 Þ jS1 j2 6¼ : Re ðS2 T2 Þ jS2 j2

ð19cÞ

4.4. Reconstruction using one image Suppose that one uses the imaging system in a single state ðS; T Þ, to form the image IOUT . Suppose, 2 further, that one can neglect the nonlinear term jTj j C in (12), leaving   o o uIN
ImðST  Þ loge IIN : IOUT ¼ IIN jS j2 þ 2Re ðST  Þ ð20Þ oX oX In Section 4.2, we required three images to solve this equation for the intensity and phase of the incident wavefield. Here, we show how one can make do with a single image if either of the following hold: (i) one knows a priori that IIN is independent of position (we term this the ‘‘pure phase case’’); or (ii) uIN is a linear function of loge IIN . An important special case of (ii) is obtained when one uses normally incident plane waves to illuminate an object composed of a single material, with the object being sufficiently thin for the projection approximation to hold. Case (i) is studied in Section 4.4.1, while case (ii) is treated in Section 4.4.2.

D. Paganin et al. / Optics Communications 234 (2004) 87–105

4.4.1. Pure phase case 0 If we assume that IIN ðX ; Y Þ ¼ IIN is independent of X and Y , then (20) becomes   o 2 0 uIN ; IOUT ¼ IIN jS j þ 2Re ðST  Þ oX which is readily solved for ouIN =oX to yield   0 IOUT =IIN
jS j2 o uIN ¼ ; Re ðST  Þ 6¼ 0: oX 2Re ðST  Þ

95

ð21Þ

ð22Þ

0 ¼ constant implies that If we keep the nonlinear term jTj j2 C in (12), then the assumption IN ðX ; Y Þ ¼ IIN  2 IOUT o o u þ jT j2 u ¼ jS j2 þ 2Re ðST  Þ : ð23Þ 0 oX IN oX IN IIN

Using the quadratic formula, this can be solved for ouIN =oX , to give qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
Re ðST  Þ  ½Re ðST  Þ2
jST j2 þ IOUT jT j2 0 IIN o u ¼ ; T ¼ 6 0: 2 oX IN jT j 2

ð24Þ

2

Now, ½Re ðST  Þ
jST j2 ¼
½ImðST  Þ , which allows us to simplify (24) to ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi   s  2 o S S IOUT u ¼
Re þ 0 ; T ¼ 6 0:  Im oX IN T T IIN jT j2

ð25Þ

We now show how to resolve the sign ambiguity in (25). The quantity under the square root sign will be strictly positive for every location on a given image. Assume that ouIN =oX is continuous, which will be the case if no quantized phase vortices are present [13]. Since (i) ouIN =oX is continuous, (ii) Re ðS=T Þ is constant and (iii) the square-root term of Eq. (25) is never zero, one concludes that (iv) either the plus or the minus sign must hold for all locations in a given image. To complete the argument, we need to show how this sign can be determined. To this end, consider a region of the image which is known a priori to have no sample in its immediate vicinity, implying that ouIN =oX ¼ 0 within this region; at any point contained within such a 0 region, Eq. (23) becomes jS j2 ¼ IOUT =IIN . Hence, when ouIN =oX ¼ 0, Eq. (25) becomes Re ðS=T Þ ¼ 1, qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 1  ½ImðS=T Þ þ jSj =jT j > 0; this implies that (i) the plus sign in (25) should be chosen if Re ðS=T Þ > 0 ; the negative sign in (25) should be chosen if Re ðS=T Þ < 0.

4.4.2. Homogeneous object case Here, we assume that uIN is a linear function of loge IIN , i.e. c uIN ¼ loge IIN þ C; ð26Þ 2 where c and C are real constants. Suppose that one uses the imaging system in a single state ðS; T Þ to form the image IOUT . Assume, further, that one can neglect the nonlinear term in the square brackets of (11), leaving (20) as the equation to be solved. If we substitute (26) into (20), we obtain  IOUT ¼

2 jS j þ s

o oX

 IIN ;

s  cRe ðST  Þ
ImðST  Þ:

ð27Þ

As pointed out at the end of Section 4.1, equations of this form can be easily solved using the Fourier derivative theorem. Specifically, use the convention defined by Eq. (9) to take the Fourier transform of

96

D. Paganin et al. / Optics Communications 234 (2004) 87–105 2

Eq. (27) with respect to X and Y , leading to the equation FIOUT ¼ ðjSj þ isnÞFIIN . This algebraic equation can then be solved for FIIN ; the inverse Fourier transform of this expression then gives ! " # 1
1 IIN ¼ F FIOUT : ð28Þ 2 jSj þ isn Eq. (28) is the desired single-image reconstruction formula for IIN . Once IIN has been so reconstructed, it can be substituted into (26) in order to determine uIN . One is thereby led to the following formula for pffiffiffiffiffiffi reconstructing the complex wavefunction wIN ¼ IIN expðiuIN Þ which was input into the imaging system ( " #)ð1þicÞ=2 ! 1
1 wIN ðX ; Y Þ ¼ expðiCÞ F : ð29Þ FIOUT ðX ; Y Þ 2 jSj þ isn Remarkably, there is no non-degeneracy condition for Eq. (29) . This equation is a non-linear integral transform mapping the measured intensity data IOUT ðX ; Y Þ to the complex wavefunction wIN ðX ; Y Þ which resulted in this data. Note that we can set C ¼ 0, since it is a global phase factor that can be set to zero by gauge freedom. We now look at an important special case of (29), which applies when one has a ‘‘homogeneous object’’. Here, we consider an object or sample to be ‘‘homogeneous’’ when it is composed of a single material. Let 0 the homogeneous object be illuminated with normally incident plane waves, of uniform intensity IIN , which
are travelling in the z direction, and let the projected thickness of material in the z direction be T ðX ; Y Þ. If the object is sufficiently thin for the projection approximation to be valid, then the intensity IIN ðX ; Y Þ and phase uIN ðX ; Y Þ of the radiation at the exit surface of the object will be well approximated by: uIN ¼
kdT
þ u0IN ;

0 IIN ¼ IIN expð
lT
Þ:

ð30Þ

Here, k ¼ 2p=k, k is the wavelength of the radiation, n ¼ 1
d is the real part of the refractive index of the material, l is its absorption coefficient, and u0IN is a constant phase factor which may be set to zero without loss of generality. The first member of (30) states that the phase shift, imparted on the radiation by the sample, is proportional to the projected thickness of the sample; the second member of (30) is the Beer– Lambert law of absorption. The ‘‘homogeneous object’’ Eqs. (30) imply that uIN and IIN obey Eq. (26), with c¼

2kd ; l

C ¼ u0IN


kd loge I0IN : l

ð31Þ

This allows us to make use of the reconstruction formulae Eqs. (28) and (29), to obtain IIN ðX ; Y Þ and wIN ðX ; Y Þ from IOUT ðX ; Y Þ. If one is interested in the projected thickness T
ðX ; Y Þ of material in the sample, then substitution of the reconstructed intensity IIN into the second member of (30) gives an equation which is trivially solved for T
ðX ; Y Þ. Making use of (31), the resulting formula is 91 0 8  = < 1 1 I ðX ; Y Þ OUT A: h iF ð32Þ T
ðX ; Y Þ ¼
loge @F
1 0 : jS j2 þ in 2kd Re ðST  Þ
ImðST  Þ ; l IIN l

Eq. (32) is a single-image algorithm for reconstructing the projected thickness T
ðX ; Y Þ of a single-material sample which has been imaged using a coherent shift-invariant linear imaging system with a linear transfer function. In general, the reconstruction depends on seven real parameters, namely the two complex pa0 rameters ðS; T Þ characterising the imaging system, the intensity IIN of the incident plane wave, the absorption coefficient l of the sample, and the real part 1
d of the refractive index of the sample. Note, 0 however, that only three real parameters enter into the reconstruction formula (32), namely l, IIN jS j2 and

D. Paganin et al. / Optics Communications 234 (2004) 87–105

97

0 IIN ½2kdl
1 Re ðST  Þ
ImðST  Þ. If desired, one can substitute the result of (32) into (30), in order to reconstruct both the intensity and phase of the complex disturbance wIN ðX ; Y Þ of the wavefield incident on the imaging system. Further, we note that if only uIN ðX ; Y Þ and IIN ðX ; Y Þ are to be reconstructed, then d=l is the only information one needs regarding the single-material sample being imaged. Next we show how, given only a single recorded image IOUT ðX ; Y Þ, one can obtain a grey-scale image of 0 T
ðX ; Y Þ as a function of a single real parameter s. While the relevant parameters S; T ; l; d; IIN must be known if one is to perform a quantitative reconstruction of T
ðX ; Y Þ, one is often only interested in a greyscale map of the same. In this context, note that linear-palette grey-scale maps of a given function are unchanged by (i) the addition of a constant to the function and (ii) multiplication of the function by a given positive number, if the images are scaled such that black corresponds to the minimum value of the function and white corresponds to its maximum value. Utilising this, we re-arrange (32) into the form 8 9 < = FI ðX ; Y Þ OUT 2 0 h
1 i lT
ðX ; Y Þ
loge ðIIN : ð33Þ jS j Þ ¼
loge F
1 : 1 þ in 2kdl Re ðST  Þ
ImðST  Þ ; jS j2

Note that a grey-scale image of the left side of (33) will be the same as a grey-scale image of T
ðX ; Y Þ itself. In this context, denote the left side of (33) by T
ðX ; Y Þ, this being the projected thickness up to irrelevant multiplicative and additive constants. Also, let all the physical parameters in square brackets be denoted by a single real parameter s: s

2kdl
1 Re ðST  Þ
ImðST  Þ 2

jS j

;

ð34Þ

which is a special case of the second member of Eq. (27) . Therefore (33) has the following compact form


1 FIOUT ðX ; Y Þ
T ðX ; Y Þ ¼
loge F : ð35Þ 1 þ ins The reconstructed grey-scale map T
ðX ; Y Þ of the projected thickness of the sample is thus a function of the single image IOUT ðX ; Y Þ, together with the real parameter s. At this point, we introduce an analogy with optical systems, such as microscopes and cameras, which possess a variable focus setting. Here, the state of a hardware imaging system is tunable via a single continuous parameter, which is varied until one judges ‘‘by eye’’ that the resulting image is in focus. Similarly, the algorithm in (35) – which takes the output IOUT ðX ; Y Þ of the imaging system as an input and produces a grey-scale image T
ðX ; Y Þ of the projected thickness as an output – depends on a single continuous parameter s; which can be adjusted ‘‘by eye’’ until the image is deemed to be ‘‘in focus’’. This value of s will be that for which the characteristic ‘‘derivative’’ contrast just disappears (cf. Eq. (27)). The parameter s in (35) may therefore be viewed as a generalised focal parameter in a digital optical system [14], i.e., one which uses a mathematical lens rather than a hardware lens to bring an image into focus. Our reconstruction of the projected thickness, which may be viewed as a decoding of the coded DEI image, allows us to generate the image which would have resulted if s were equal to zero in the imaging system. To take this point further, note that once the complex wavefunction incident on the imaging system has been reconstructed using Eq. (29), one can emulate in software the action of any subsequent imaging system for which an associated mathematical transform can be written down [15]. An attractive feature of (35) is that it can be used to provide a grey-scale reconstruction of T
ðX ; Y Þ, given a single image obtained using a coherent shift-invariant linear imaging system with a linear transfer function; no other knowledge about the state of the system is necessary, given that the single real parameter s in the reconstruction can be deduced using the method described above. This numerically simple means of reconstruction may therefore be applied to images taken in real time; this will permit, for example, the

98

D. Paganin et al. / Optics Communications 234 (2004) 87–105

analysis of movies of the evolution of single-material objects taken using both DEI (see Section 5) and DIC imaging (see Section 6). Lastly, note that the single-material assumption can be extended slightly, with no change in the mathematics. If, rather than a single-material sample, one has a ‘‘quasi-homogeneous’’ sample consisting of a single-material substrate of approximately constant (or slowly varying) projected thickness along the direction of the optic axis, into which is embedded a different material (see Fig. 1), then all of the reconstruction formulae of this section hold true. One need only replace the absorption coefficient l with the difference between the absorption coefficients for sample and substrate (‘‘effective absorption coefficient’’), and replace the real part of the refractive index 1
d with the difference between the real parts of the refractive indices for sample and substrate (‘‘effective refractive index’’). Use will be made of this fact in Sections 5.2 and 6.

5. Application: diffraction enhanced imaging using hard X-rays The work of Lang [9] contains the central idea behind the technique of diffraction enhanced imaging. His ‘‘projection topograph’’, which is formed by imaging the reflection of a well-collimated X-ray beam from an imperfect crystal, is extremely sensitive to dislocations and other imperfections in the crystal. In terms of geometric optics, this sensitivity may be roughly pictured as due to the fact that perturbations introduced by the crystalÕs defects strongly affect the intensity of a reflected ray. One may view LangÕs imperfect crystal as consisting of the sum of a perfect crystal (‘‘analyzer’’) and its deviations from perfection (‘‘sample’’). This leads directly to the idea of separating the sample from the analyzer, which is the essence of the idea published by F€ orster and colleagues in 1980 [16]. Here, appropriately filtered radiation passes through a sample of interest, with the resulting perturbed wavefield being diffracted from a perfect crystal – see Fig. 1. The exquisite sensitivity, of the crystalÕs local reflection coefficient to transverse phase gradients in the beam incident upon it, results in strong phase contrast in the registered image. This elegant idea has been further developed by a number of workers [8,17–20]. We adopt the widely used terminology introduced by Chapman et al. [8], and refer to such crystal-based imaging as ‘‘diffraction enhanced imaging’’ (DEI). Here, we consider the inverse problem of DEI [8,21–23] by applying the formalism of Section 4 to diffraction enhanced imaging. Section 5.1 outlines the theory of quantitative phase-amplitude reconstruction using one or more diffraction enhanced images, with particular emphasis being paid to (i) obtaining additional correction terms to the geometric-optics Chapman equations of diffraction-enhanced imaging [8] by taking the Takagi wave-optics equations [24] as a starting point, and (ii) applying the single-image reconstruction algorithm (35) to a diffraction enhanced image of a quasi-homogeneous object. Section 5.2 gives an experimental demonstration of these ideas, while Section 5.3 is devoted to a stability analysis of the single-image reconstruction algorithm. 5.1. Theory With reference to Fig. 1, consider normally incident plane X-rays that are weakly deformed by passage through a given sample, and which then strike the inclined surface of a perfect crystal. The crystal is oriented such that, in the absence of an object, the plane wave strikes the crystal at or sufficiently near to the Bragg angle for a given reflection, such that the plane wave is reflected from the crystal. With the sample present, the intensity of the reflected wave is recorded using a two-dimensional imaging device. Both the complex scalar wavefield incident on the crystal, and the wavefield reflected from the crystal, can be represented as deformed plane waves. These deformed plane waves can be written as a product of two factors, namely: (i) a rapidly varying complex amplitude representing the plane wave which would exist in the absence of the object, and (ii) slowly varying envelopes which quantify the perturbation induced on the

D. Paganin et al. / Optics Communications 234 (2004) 87–105

99

plane waves by their passage through the sample, where ðx; yÞ denote Cartesian coordinates over the surface of the crystal. This is an example of a coherent shift-invariant linear imaging system, as described in Section 3. If the refraction angles (of the coherently scattered X-rays passing through the object) are sufficiently small, then one may safely linearize the transfer function associated with this imaging system, as described by Eq. (5). ^ A in Eq. (1), namely the transfer function, which For this system, the Fourier transform of the operator O relates the slowly varying component wIN ðx; yÞ of the wavefield incident on the crystal to the slowly varying component wOUT ðx; yÞ reflected from the crystal, is known as the Takagi propagator [10]. Note that ^ A , for the just-introduced im‘‘propagator’’ is a synonym for ‘‘transfer function’’. Denote the operator O ^ TAK . Transform to a rotated coordinate system (X ; Y ) such that the Fourier transform of aging system, by O ^ TAK is independent of the Fourier space coordinate g (note that this coordinate is dual to Y ; see Fig. 2). O Evidently, the formalism of Section 3 may be applied to this system, without needing to invoke the specific functional form of the Takagi propagator for a given crystal. Note that, in this application, the direction X is the coordinate along the crystal surface which lies in the diffraction plane; n is dual to X . The direction Y is along the crystal surface but perpendicular to the diffraction plane; g is dual to Y . ^ TAK . This is related ~ TAK ðn; gÞ ¼ K ~ TAK ðnÞ  FO Bearing Eq. (4) in mind, write the Takagi propagator as K to the rocking curve RðnÞ of the crystal by ~ TAK ðnÞj2 ; Rðn0 þ nÞ ¼ jK

ð36Þ

where n0 is a parameter related to both the crystal geometry and its inclination with respect to the incident beam (see [10], where this parameter is called ‘‘A’’). We now linearize the Takagi propagator about the point ðn; gÞ ¼ ð0; 0Þ in Fourier space, which corresponds to An ¼ Ag ¼ 0 (see Fig. 2). Next, we need to calculate the ordered pair ðS; T Þ ¼ ðSTAK ; TTAK Þ which completely characterizes the imaging system. Evidently,   o ~  : TTAK  K ðnÞ ð37Þ TAK  on n¼0

Substitute this into the second member of (8) in order to construct STAK : ~ TAK ðn ¼ 0Þ; STAK  K

ð38Þ

where we have used the fact that An ¼ 0. Now that we have the ordered pair ðS; T Þ ¼ ðSTAK ; TTAK Þ, we see that Eq. (11) becomes: !   o   ðSTAK TTAK Þ oXo uIN
Im SoTAK TTAK loge IIN jSTAK j2 þ 2Re oX n  2  2 IDEI ¼ IIN ; ð39Þ þjTTAK j2 oXo uIN þ 14 oXo loge IIN where IOUT ðX ; Y Þ has now been written as the DEI image IDEI ðX ; Y Þ, and the complex numbers ðS; T Þ ¼ ðSTAK ; TTAK Þ are as given in Eqs. (37) and (38). For notational clarity, we have silently assumed that the same coordinate system (X ; Y ) is used to measure displacements in the plane of the crystal and in the plane of the imaging detector; this assumption is readily lifted, if required (e.g. for asymmetric reflections). Substitute Eqs. (37) and (38) into (39), then make use of both (36) and o=on of Eq. (36) evaluated at ðn; gÞ ¼ ð0; 0Þ, to arrive at:  3  o Rðn0 Þ þ ono Rðn þ n0 Þ u oX IN 7 6  ðn;gÞ¼ð0;0Þ    7 6  7 6 o ~ o ~ K
Im K ð0Þ ðnÞ log I  7: 6 TAK TAK IN e ¼ IIN 6 on oX 7 n¼0 6    2 n o7 5 4 o    2 2  o ~ TAK ðnÞ  þ on K u þ 14 oXo loge IIN  oX IN 2

IDEI

n¼0

ð40Þ

100

D. Paganin et al. / Optics Communications 234 (2004) 87–105

This equation is a special case of the general expression (11). One may therefore use any of the methods outlined in Section 4, in order to quantitatively reconstruct the intensity IIN and phase uIN of the slowly varying component of the wavefield striking the crystal, given as data the DEI images obtained using one or more orientations of the analyzer crystal. Without further approximation, the methods of Section 4.1 can be used to reconstruct the ‘‘apparent absorption image’’ IIN , the ‘‘apparent refraction image’’ ouIN =oX , and the phase map uIN , given four images of the incident wavefield, obtained using four orientations of the analyzer crystal. If we are prepared to discard the non-linear term    ~ TAK ðnÞ=on j2 ðouIN =oX Þ2 þ 1 ðo loge IIN =oX Þ2 joK 4 n¼0 from (40), then the methods of Section 4.2 allow us to perform the same reconstructions using three images. If we approximate further by keeping only the first two terms in the square brackets of (40), then we arrive at the Chapman equations [8]: " #   o o  Rðn þ n0 Þ uIN : IDEI ¼ IIN Rðn0 Þ þ ð41Þ on n¼0 oX Eq. (19) of Section 4.3 then reduces to those which have been published previously [8], if one assumes that: (i) the rocking curve is symmetric, and (ii) the two analyzer crystal positions are symmetrically placed about the peak of the rocking curve. For a quasi-homogeneous object (as defined at the end of Section 4), the discarding of Eq. (40)Õs non-linear term implies that the reconstruction may be performed using a single image; in particular, the analysis of Section 4.4.2 can be used to obtain a reconstructed grey-scale map T
ðX ; Y Þ of the projected thickness of the quasi-homogeneous sample, without any detailed knowledge of the structure or orientation of the analyzer crystal. This single-image reconstruction takes into account a term linear in the derivative of the logarithm of intensity, which is missing from (41). Finally, we note that when either or both of the transverse absorption gradient and the transverse intensity gradient are sufficiently strong (e.g., for X-ray wavefields passing through the edges of sharp objects), then the additional terms in (40) may make this equation more suitable for quantitative DEI analyses than Eq. (41). 5.2. Experiment DEI images were obtained at the SYRMEP beamline of the Elettra Synchrotron Light Laboratory in Trieste, Italy. Monochromatic radiation at 20 keV was selected using a Si(1 1 1) channel cut monochromator located 12 m upstream of the imaging system. The distance from the sample to the analyzer crystal was approximately 70 cm, with the distance from the crystal to the detector at around 75 cm. Images were obtained by scanning an image plate vertically through the beam, with the image plate subsequently being read out on a Fuji BAS 1800 at 50 lm resolution. See [25] for a more detailed description of the experimental setup. As a well-characterised quasi-homogeneous test sample, use was made of a TOR(MAM)Ò mammographic test phantom (Department of Medical Physics, Leeds University, UK). This phantom contains features which include 10 mm long fibres (nylon fishing line), together with 3 mm diameter plastic disks of varying densities. The phantom, which is just over 1 cm thick, had a 3.3 cm thick slab of polymethylmethacrylate (PMMA) placed against it, so as to more realistically emulate the imaging of a full thickness of breast tissue [25]. Using the previously described setup, the DEI images shown in Figs. 3(a) and (b) were recorded. Fig. 3(a) shows a 350  350 pixel DEI sub-image of a geometric arrangement of four nylon fibres, each of which is approximately 10 mm long. Fig. 3(b) shows a 300  300 pixel DEI sub-image of three plastic disks,

D. Paganin et al. / Optics Communications 234 (2004) 87–105

101

Fig. 3. (a) DEI sub-image of four nylon fibre inclusions in a TOR(MAM)Ò breast phantom. Image size: 350  350 pixels. (b) DEI subimage of three disks, two of which are denser than the third, in the TOR(MAM)Ò breast phantom. Image size: 300  300 pixels. (c) Grey-scale reconstruction of the projected thickness corresponding to the image in (a), obtained using Eq. (35), using (a) as the only input data. (d) Grey-scale reconstruction of the projected thickness corresponding to the image in (b), obtained using Eq. (35). The scale bar in all images is 2 mm in length.

each of which have 3 mm diameter, and two of which are much denser than the remaining disk. Note that the experimentally obtained images have been rotated by 90°. The two mechanisms of contrast formation, i.e., absorption contrast and refraction contrast, are clearly evident in Fig. 3(a) . The weak contrast of the horizontal fibre is, with the exception of the contrast at its ends, due exclusively to absorption. This absorption contrast is present for all fibres, but for non-horizontal fibres the additional refraction (‘‘derivative’’) contrast provided by the analyzer crystal is dominant. Similarly, the ‘‘derivative’’ contrast in the image of the disks is far stronger than the barely discernible absorption contrast due to the bulk of the disk. The mathematical origin of the derivative contrast is evident in the derivation leading to Eq. (40), while its physical origin lies in the extreme angle-specificity of the crystalÕs rocking curve RðnÞ [8,25]. A short computer program was written to implement the single-image reconstruction algorithm in Eq. (35), this being applied to the images in Figs. 3(a) and (b). Making use of the fast Fourier transform in its

102

D. Paganin et al. / Optics Communications 234 (2004) 87–105

implementation, each reconstruction took approximately 2 seconds on a Pentium III 1 GHz desktop computer running the Interactive Data Language (IDL) version 5.6. We then varied the single parameter s characterising this virtual imaging system, in analogy with the ‘‘focus knob’’ on a microscope, until the characteristic black-white contrast of the DEI images just disappeared (see Section 4.4.2). This led to the grey-scale maps of the projected thickness T
ðX ; Y Þ (up to irrelevant multiplicative and additive constants) which are shown in Figs. 3(c) and (d). Note that no attempts were made to stabilise the reconstruction by well-known techniques such as regularisation [26,27], pre- and post-smoothing of the data, solutions of minimum norm [27], etc. Our implementation of the single-image DEI reconstruction algorithm in Eq. (35), which makes simultaneous use of information contained in both absorbed and refracted X-rays, yields quantitatively correct grey-scale maps which are a direct representation of the projected thickness of the sample, unlike the raw DEI images themselves. In particular, the reconstruction removes the imaging systemÕs bias for sample edges and variations along a single direction, both of which serve to mask bulk structure (such as the interior of the disks) if that structure is slowly varying along a direction orthogonal to the diffraction plane. The algorithm both makes use of, and serves to redress, the DEI imageÕs strong bias for edge information and structures orthogonal to the diffraction plane. Since both refracted and absorbed X-rays are used in the reconstruction, one is able to make do with smaller doses compared to those which would be needed if absorption were the sole contrast mechanism. Indeed, the reconstructed images of the sampleÕs projected thickness can be used to calculate the absorption radiographs one would have obtained if (i) the dose were higher and (ii) no wideangle scattered radiation was present (this latter point is particularly important in light of the large degree of multiple scatter from thick samples). Further work is clearly needed to thoroughly investigate the possible utility, for improved methods of tumour diagnosis, of this and other methods for quantitative DEI-image analysis (cf. [25]). Simplicity is an obvious advantage of our single-image DEI reconstruction technique: (i) no provision need be made for the taking of multiple DEI images and their subsequent alignment; (ii) no experimental parameters are required in the reconstruction, which takes a single DEI image as the only input. Indeed, the single-image method is ideally suited to the quantitative analysis of real-time DEI or DIC movies, and as such may be of utility in contexts like the quantitative real-time study of plasmas using X-ray DEI or DIC imaging, and the quantitative visible-light DIC imaging of biological cells and cold atomic clouds. 5.3. Stability analysis Despite its stability with respect to a significant amount of noise in the input data, there is a clear streaking artefact in the reconstructions of Figs. 3(c) and (d) . The origin of this artefact can be understood from both a physical and a mathematical point of view, as we now show. (a) The refraction contrast provided by the analyzer crystal is sensitive to ray deviations in one direction only, namely the direction parallel to the analyser diffraction plane (i.e., the X direction). Therefore, in passing from a given sample to its DEI image, X -varying spatial information is greatly amplified with respect to Y-varying spatial information, where Y is the image-plane axis perpendicular to X . Thus, in the inverse problem of going from the DEI image to the projected thickness of the sample, X -varying spatial information must be greatly suppressed with respect to Y-varying spatial information. Thus the noise-induced deviations in the data introduce a ‘‘streak’’ artefact into the reconstruction, which is much more strongly varying in the Y than the X direction. (b) The single-image reconstruction algorithm (35) displays maximum sensitivity to noise-induced perturbations in the data, for spatial frequencies expðirgÞ; r 2 R, corresponding to the line n ¼ 0 in Fourier space, for it is along this line that 1=ð1 þ insÞ is greatest in modulus. Since such spatial frequencies have no variation in the X direction, the resulting artefact is dominated by the observed streaking along this

D. Paganin et al. / Optics Communications 234 (2004) 87–105

103

direction. Significantly, since 1=ð1 þ insÞ in (35) is well behaved for any coordinate ðn; gÞ in Fourier space, our single-image phase retrieval algorithm is stable.

6. Application: quantitative differential interference contrast microscopy (DIC) Here, we show how the ideas of this paper may be used to render quantitative the images which are taken using a differential interference contrast microscope, provided that the (quasi) homogeneous object to be imaged is sufficiently weak, in a sense which shall soon be made precise. An ability to quantitatively analyse visible-light DIC microscope images is of utility in the context of measuring the dry mass of biological objects in vivo [28], and in the quantitative analysis of real-time biological systems imaged, without staining, through a differential interference contrast microscope [29]. These ideas may also be useful in the context of X-ray DIC microscopy [30,31]. Consider plane monochromatic waves travelling through an object of projected thickness T
ðx; yÞ, which is embedded in a substrate of some form. In the context of optical microscopy, the substrate will typically consist of the slab of fluid between the upper surface of the microscopeÕs glass slide and the underside of the coverslip on top of this fluid. If the object in the fluid is sufficiently weakly scattering, and the substrate is of uniform thickness, then the wavefunction at the exit surface of the sample may be calculated using the projection approximation (cf. Eq. (30)), to yield:    qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 wðx; yÞ ¼ expð
DlT
ðx; yÞÞ expð
ikDdT
ðx; yÞÞ ¼ exp
T
ðx; yÞ Dl þ ikDd : ð42Þ 2 Here, Dl is the difference between the absorption coefficient lo of the object and the absorption coefficient ls of the substrate, Dd is the difference between the real part ns ¼ 1
ds of the refractive index of the substrate and the real part no ¼ 1
do of the refractive index of the object, and we have discarded an irrelevant multiplicative complex constant. The differential interference contrast image is formed by taking the exit-wavefunction in (42), and then interfering it with a slightly displaced copy of itself that has had a constant phase shift /0 imposed upon it. Assuming this slight displacement to be by an amount dx in the x direction, the resulting ‘‘DIC’’ wavefunction wDIC is: wDIC ðx; yÞ ¼ wðx; yÞ þ expði/0 Þwðx
dx; yÞ       1 1 Dl þ ikDd þ expði/0 Þ exp
T
ðx
dx; yÞ Dl þ ikDd : ¼ exp
T
ðx; yÞ 2 2

ð43Þ

Assume the displacement dx to be sufficiently small to permit the first-order Taylor approximation T
ðx
dx; yÞ  T
ðx; yÞ
dxoT
ðx; yÞ=ox, which may then be substituted into (43). If we take the square modulus of the resulting expression, we arrive at the following form for the intensity IDIC of the DIC image: i h i
IDIC ðx; yÞ ¼ exp
Dl T ðx; yÞ 1 þ exp Dldx oT
ðx; yÞ=ox   h i

1 Dl dx oT
ðx; yÞ=ox cos k Dd dx oT
ðx; yÞ=ox þ /0 : þ 2 exp 2 h

ð44Þ

Expand the cosine term using the relevant double-angle formula, and then assume the projected thickness T
ðx; yÞ to be sufficiently weak that one may safely replace all transcendental functions of the projected thickness, and its partial derivative with respect to x, with first-order Taylor expansions. In the resulting expression, keep only those terms which are linear in T
ðx; yÞ and its partial derivative with respect to x, leaving:

104

D. Paganin et al. / Optics Communications 234 (2004) 87–105

1 a
IDIC ðx; yÞ ¼ ½Dla þ bo=oxT
ðx; yÞ; 2

a  1 þ cos /0 ;

1 b  kDddx sin /0
Dldxa: 2

ð45Þ

This equation is mathematically identical in form to Eq. (27) of the main text, and may therefore be solved for T
ðx; yÞ using the same method. As a corollary to the mathematical identity just mentioned, we note that, for the weak T
case to which the present calculation restricts itself, DIC and DEI images are indistinguishable from one another, despite the difference in the physical mechanisms of image formation. For a concrete demonstration of this fact, one may compare the DEI image in Fig. 3(a) of the present paper to the DIC image in Fig. 6 of reference [29].

7. Conclusion We have considered the problem of quantitative phase-amplitude retrieval from images obtained using a coherent shift-invariant linear imaging system with a linear transfer function. Special cases of such systems include those used for differential interference contrast microscopy with visible light or X-rays, and diffraction enhanced imaging with hard X-rays. We gave several algorithms for reconstructing both the phase and amplitude of a given two-dimensional coherent wavefield, given as input data one or more images of such a wavefield which may be formed by different states of the coherent linear imaging system. We found that, when an object to be imaged consists of a single material embedded within a substrate of approximately constant projected thickness, the phase-amplitude reconstruction can be performed using a single image. This single-image algorithm was successfully applied to the problem of quantitatively imaging inclusions in a breast tissue phantom, using a single synchrotron DEI image of the same. In particular, the reconstructed grey-scale images properly represent the bulk of the sample, in contrast to DEI images which greatly emphasise sharp edges. The ability to obtain reconstructions, from a single DEI image, is a particularly attractive feature of the method which will allow quantitative DEI reconstructions to be performed on real-time data. Lastly, correction terms were obtained to the Chapman equations of diffraction enhanced imaging: these correction terms take into account both transverse absorptive gradients, and the non-linear contribution of both transverse phase gradients and transverse absorptive gradients.

Acknowledgements David Paganin and Konstantin Pavlov acknowledge support from the Australian Research Council. Marcus Kitchen acknowledges the receipt of an Australian Postgraduate Award. The authors acknowledge stimulating discussions with Michael Morgan, Karen Siu and Andrew Pogany. The experimental data, used in the reconstructions given in this paper, was obtained with the kind assistance of Chris Hall, Alan Hufton, and the superb staff of the SYRMEP beamline at Elettra.

References [1] [2] [3] [4] [5] [6]

M. Born, E. Wolf, Principles of Optics, sixth ed., Pergamon Press, Oxford, 1993. F. Zernike, Physica 9 (1942) 686. M. De Graef, Magnetic Imaging and its Application to Materials: Experimental methods in the Physical Sciences vol. 36 (2001) 27. D. Gabor, Nature 161 (1948) 777. G. Nomarski, A.R. Weill, Rev.Metall. 2 (1955) 121. J.R. Meyer-Arendt (ed.), Selected papers on Schlieren optics, SPIE Milestone series vol. MS61, The International Society for Optical Engineering, Bellingham, Washington, 1992. [7] W. Coene, G. Janssen, M. Op de Beeck, D. Van Dyck, Phys. Rev. Lett. 69 (1992) 3743.

D. Paganin et al. / Optics Communications 234 (2004) 87–105

105

[8] D. Chapman, W. Thomlinson, R.E. Johnston, D. Washburn, E. Pisano, N. Gm€ ur, Z. Zhong, R. Menk, F. Arfelli, D. Sayers, Phys. Med. Biol. 42 (1997) 2015. [9] A.R. Lang, Acta Crystallogr. 12 (1959) 249. [10] T.E. Gureyev, S.W. Wilkins, Il Nuovo Cimento 19D (1997) 545. [11] M. Nieto-Vesperinas, Scattering and Diffraction in Physical Optics, Wiley, New York, 1991. [12] A. Papoulis, Systems and Transforms with Applications in Optics, McGraw-Hill, New York, 1968. [13] L. Pismen, Vortices in Nonlinear Fields: From Liquid Crystals to Superfluids, from Non-equilibrium Patterns to Cosmic Strings, Oxford University Press, Oxford, 1999. [14] L. Yaroslavsky, M. Eden, Fundamentals of digital optics: Digital signal processing in optics and holography, Birkh€auser, Boston, 1996. [15] D. Paganin, T.E. Gureyev, S.C. Mayo, A.W. Stevenson, Ya.I. Nesterets, S.W. Wilkins, X-ray Omni Microscopy, J. Microsc., 2004 (in press). [16] E. F€ orster, K. Goetz, P. Zaumseil, Krist. Tech. 15 (1980) 937. [17] V.A. Somenkov, A.K. Tkalich, S. ShilÕstein, Sov. Phys. Tech. Phys. 3 (1991) 1309. [18] V.N. Ingal, E.A. Beliaevskaya, J. Phys. D 28 (1995) 2314. [19] T.J. Davis, D. Gao, T.E. Gureyev, A.W. Stevenson, S.W. Wilkins, Nature 373 (1995) 595. [20] T.J. Davis, T.E. Gureyev, D. Gao, A.W. Stevenson, S.W. Wilkins, Phys. Rev. Lett. 74 (1995) 3173. [21] E. Pagot, P. Cloetens, S. Fiedler, A. Bravin, P. Coan, J. Baruchel, J. H€artwig, Appl. Phys. Lett. 82 (2003) 3421. [22] V.A. Bushuev, A.A. Sergeev, Tech. Phys. Lett. 25 (1999) 83. [23] C. Preza, J. Opt. Soc. Am. A 17 (2000) 415. [24] S. Takagi, J. Phys. Soc. Jpn. 26 (1969) 1239. [25] R.A. Lewis, K.D. Rogers, C.J. Hall, A.P. Hufton, S. Evans, R.H. Menk, G. Tromba, F. Arfelli, L. Rigon, A. Olivo, A. Evans, S.E. Pinder, E. Jacobs, I.O. Ellis, D.R. Dance, in: L.E. Antonuk, M.J. Yaffe (Eds.), Medical Imaging 2002: Physics of Medical Imaging, Proc. SPIE vol. 4682 (2002) 286. [26] A.N. Tikhonov, Dok. Akad. Nauk. SSSR 3 (1963) 501. [27] R. Kress, Linear Integral Equations, Springer, Berlin, 1989. [28] E.B. van Munster, E.K. Winter, J.A. Aten, J. Microsc. 191 (1998) 170. [29] E.B. van Munster, L.J. van Vliet, J.A. Aten, J. Microsc. 188 (1997) 149. [30] T. Wilhein, B. Kaulich, E. Di Fabrizio, F. Romanato, S. Cabrini, J. Susini, Appl. Phys. Lett. 78 (2001) 2082. [31] B. Kaulich, T. Wilhein, E. Di Fabrizio, F. Romanato, M. Altissimo, S. Cabrini, B. Fayard, J. Susini, J. Opt. Soc. Am. A 19 (2002) 797.