WO2003034010A1

WO2003034010A1 - Phase determination of a radiation wavefield

Info

Publication number: WO2003034010A1
Application number: PCT/AU2002/001398
Authority: WO
Inventors: Brendan E Allman; Keith Nugent
Original assignee: Iatia Imaging Pty Ltd
Priority date: 2001-10-16
Filing date: 2002-10-15
Publication date: 2003-04-24
Also published as: AUPR830801A0

Abstract

The phase of radiation wavefield is retrieved by solving the transport of intensity equation. The rate of change of intensity, orthogonal to a surface extending across the wavefield, is first determined (eg, by measuring intensities at two separated planes). This rate of change is subjected to the computational process of taking an integral transform, multiplying by a filter corresponding to the inversion of a differential operator, and taking the inverse integral transform. The result is multiplied by a function of the intensity over the surface, and subjected to a similar computational process, to obtain a measure of the phase over the surface. The filters have a form based on details of the optical system used to acquire the intensity data, such as numerical apertures and spatial frequencies.

Description

PHASE DETERMINATION OF A RADIATION WAVEFIELD

Field of the Invention

This invention relates to the determination of phase of a radiation wavefield to enable phase images of objects to be produced. This invention is an improvement to the method for determining the phase of a radiation wavefield disclosed in International Patent Application No. PCT/AU99/00949. The contents of this International application are incorporated into this specification by this reference. As in the above International application, the term radiation wavefield is intended to include all forms of radiation that propagate in a wavelike manner, including but not limited to x-rays, visible light and electrons.

Background of the Invention

The method disclosed in the above International application enables phase data relating to the radiation wavefield to be determined from which a phase image can be constructed. The determination of the phase data relating to the radiation wavefield also enables other image modalities to be produced such as DIC, Zernike, Hoffman Contrast Images and Darkfield images. Phase images are an important tool to microscopists because they enable detail of some objects to be ascertained which are not available in conventional contrast images. In particular, fine edge detail in transparent structures such as cells and other biological samples may be more visible when phase images are taken rather than conventional contrast images.

Furthermore, the phase data, as described in Australian Provisional Application No. PR5928, also enables the other modalities referred to above to be determined in software rather than by optics in a microscope.

The method for determining the phase of a radiation wavefield according to the above International application requires the so-called transport of intensity equation to be solved. This is achieved by producing a representative measure of the rate of change of intensity of the radiation wavefield over a selected surface extending generally across the wavefield, producing a representative measure of intensity of the radiation wavefield over the selected surface, transforming the measure of rate of change to produce a first integral transform representation and applying to that representation a first filter corresponding to the inversion of a first differential operator. The inverse of the first integral transform is then applied to the first modified integral transform representation and a correction based on the measure of intensity over the selected surface to the untransformed representation is then applied. The correct untransformed representation is then transformed to produce a second integral transform representation and a second filter corresponding to the inversion of a second differential operator is applied to that representation. An inverse of the second integral transform is then applied to the second modified integral transform representation to produce a measure of phase of the radiation wavefield across a selected plane.

The form of the first and second differential operators according to the teachings of the above International application has the form

in real space and incorporates a term with the form

k² + ²

in its Fourier representation. In the above equation, k is the spatial frequency and α is an arbitrary constant which is included so that when k approaches 0, the operator does not diverge because of the inclusion of the factor ². The inversion of a first differential operator which is frequency dependent will, if the factor α is not included, tend to highlight or exaggerate low frequencies relative to high frequencies thereby masking the high frequencies which will cause a degradation in sharpness of the image obtained from the phase data over what may have otherwise been achieved.

Since fine edge detail in transparent samples is likely to be high frequency in nature, the minimisation of the high frequencies therefore results in a loss of information in the image produced from the phase data and therefore lack of clarity in the compiled image.

In order to prevent the exaggerating of the low frequencies relative to the high frequencies, the factor α in the above equation decreases the value of the above operator when frequency is small so that the low frequency signals do not swamp or become over-exaggerated compared to the high frequency signals which are likely to contain most of the information of interest.

Summary of the Invention

The objection of the present invention is to provide further improvements to the method and, in particular, to the form of the differential operator used in the method.

The invention may be said to reside in a method of quantitative determination of a phase of a radiation wavefield including the steps of:

(a) producing a representative measure of the rate of change of intensity of said radiation wave field over a selected surface extending generally across the wave field;

(b) producing a representative measure of intensity of said radiation wave field over said selected surface;

(c) transforming said measure of rate of change of intensity to produce a first integral transform representation and applying to said first integral transform representation a first filter corresponding to the inversion of a first differential operator reflected in said measure of rate of change of intensity to produce a first modified integral transform representation; (d) applying an inverse of said first integral transform to said first modified integral transform representation to produce an untransformed representation;

(e) applying a correction based on said measure of intensity over said selected surface to said untransformed representation;

(f) transforming the corrected untransformed representation to produce a second integral transform representation and applying to said second integral transform representation a second filter corresponding to the inversion of a second differential operator reflected in the corrected untransformed representation to produce a second modified integral transform representation;

(g) applying an inverse of said second integral transform to said second modified integral transform representation to produce a measure of phase of said radiation wave field across said selected plane; and

(h) wherein at least one of the first or second differential operator has a form based on an optical system used to acquire the radiation for producing the representative measure of the rate of change of intensity of the radiation wave ield over the selected surface extending generally across the wavefield.

Because the form of the differential operator is based on the actual optical system, the form of the differential operator is improved and therefore provides better results because the operator is based on the optical system. It should be understood that the reference to optical system in this specification includes conventional glass or plastic lenses and like imaging elements, as well as magnetic fields or electric fields which are used to condition radiation in the form of electrons in, for example, electron microscopes.

Preferably both the first and second differential operators have a form based on the optical system.

Preferably the first and second integral transforms are produced using a Fourier transform.

In one embodiment of the invention the differential operators have the form:

V ^"

T_P + ² where,

T_P(p) = 2πiόzjηT?^{' )}(p,η)dη

and

and wherein p is radial lateral spatial frequency, the η longitudinal spatial frequency and όz is the defocus distance in the plane of the object. Also

ε NA condenser

NA objective where NA_condensor and NA_obJectiw are respectively the numerical aperture of the condensor and the objective (These are settings and dimensions on the microscope) . p_obj is the maximum spatial frequency accepted by the objective.

The invention may also be said to reside in a computer program for quantitative determination of the phase of a radiation wavefield including code to perform the method steps described above.

The invention may also be said to reside in an apparatus for phase amplitude imaging of an object including: a radiation wavefield source to irradiate the object; an imaging system to focus radiation from the object to an imaging surface extending across the wavefield propagating from the object; means to produce a representative measure of radiation intensity over the imaging surface; and processing means to:

(i) produce a representative measure of intensity and a representative measure of rate of change of intensity in the direction of radiation propagation over a selected surface extending across the wave field from representative measures of radiation intensity over said image surface at said first focus and said second focus;

(ii) transform said measure of rate of change of intensity to produce a first integral transform representation;

(iii) apply to said first integral transform representation a first filter corresponding to the inversion of a first differential operator reflected in said measure of rate of change of intensity to produce a first modified integral transform representation;

(iv) apply an inverse of said first integral transform to said first modified integral transform representation to produce an untransformed representation;

(v) apply a correction based on said measure of intensity over said selected surface to said untransformed representation; (vi) transform the corrected untransformed representation to produce a second integral transform representation;

(vii) apply to said second integral transform representation a second filter corresponding to the inversion of a second differential operator reflected in the corrected untransformed representation to produce a second modified integral transform representation; and

(viii) apply an inverse of said second integral transform to said second modified integral transform representation to produce a measure of phase of said radiation wave field across said selected plane; and

(ix) wherein at least one of the first or second differential operators has a form based on the imaging system to focus the radiation from the object to the imaging surface.

Preferably the differential operator, in Fourier space, is

where,

and

Preferably both the first and second differential operators have the above form.

Brief Description of the Drawings A preferred embodiment of the invention will be described, by way of example, with reference to the accompanying drawings in which:

Figure 1A is a schematic illustration of an arrangement for determination of phase where an object is illuminated with plane wave radiation;

Figure IB is an illustration similar to Figure 1A, but with the object illuminated with point source radiation;

Figure 2 is a flow chart showing an implementation of the method of phase determination in accordance with an embodiment of this invention;

Figure 3 is a schematic illustration of an arrangement for phase amplitude microscopy using the method of the preferred embodiment of the invention; Figure 4 is a schematic drawing of an exemplary system for quantitative phase amplitude microscopy according to the preferred embodiment of the invention; and

Figure 5 is a flow chart according to the preferred embodiment.

Figures 1(a) and (b) show a schematic arrangement for phase determination in accordance with this invention where an object is illuminated by plane-wave radiation 2 or point source radiation 2 to produce transmitted beams 3.

At each point in space, an optical beam possesses two properties: intensity and pJiase. Intensity is a measure of the amount of energy flowing through each point, while phase gives a measure of the direction of the energy flow. Intensity may be measured directly, for example by recording an image on film. Phase is typically measured using interference with a "reference beam". In contrast the present method gives a non-inter erometric method for measuring phase.

Intensity can be measured over two parallel planes A, B extending across the direction of propagation of the wave field on the side remote from the incident radiation.

The present invention determines phase by providing a solution to the transport-of-intensity equation:

(1) v .(iv )= -k^-

where I is the intensity in the plane, the gradient operator in the plane is denoted V_x , k is the wave number of the radiation, and dl/dz is the intensity derivative or rate of change of intensity. Note that dl/dz is estimated from the difference of the measurements in the planes A & B shown in Figure 1, while the intensity I is given by the average of the measurements.

In order to obtain a solution to equation 1 the function A is first defined as:

(2 ) V A = N_±φ .

Thus equation (1) becomes :

Making use of a standard identity V_JL * V₁ = V_± , this may be written:

(4 ) __L ²A = -Jfea 7

where V_j_ denotes the two-dimensional Laplacian acting over the surface of the image. This equation has the following symbolic solution:

(5) A = -kVJ²d .

If the gradient operator V_x is applied to both sides of this equation, it becomes:

(6) V_xA = -JkV_iV_i ^"9_z7.

The defining equation (2) for the function A allows (6) to be transformed into:

(7) N₁ = -kV VJ²d_zI .

Dividing both sides by I then yields:

(8) V_J = -fc/-^IV_JV_±-²3_z/.

Taking the two dimensional divergence V_± • of both sides of (8), and again making use of the standard identity

V_±*V₁=V_i , then (8) becomes:

This equation has the following symbolic solution:

In order to implement a practical solution to equation (10), the following formulae are required. A suitably- well-behaved function f (x,y) may be written in the form of a two-dimensional Fourier integral:

(11) f(x, y) = \ I f(k_x,k_y)e^i(k^^y)dk_xdk_y.

A

The function f(k_x,k_y) is called the "Fourier transform" of f(x,y).

The x derivative of (11) yields:

(12) f(χ, y) = ] [ ik_χ f{ k_x,k_y ) ]e^l (k_xx+k_yy) dk dk

Hence the Fourier transform of — f (x,y) is equal to the dx

Fourier transform of f(x,y) multiplied by ik_x . Stated

differently, — = z^'F^'^F , where F denotes Fourier dx transformation and F^""1 denotes inverse Fourier transformation. Similar considerations apply to d —-f (x,y) .

By d² d² If the Laplacian V^ =—-H of (11) is obtained and dx² dy similar reasoning applied, it follows that V -_a2² _=—π F-^~ k_r F, where k² = k² +k². Thus:

Here, F denotes Fourier transformation, F^"1 denotes inverse Fourier transformation, (k_x k_y) are the Fourier variables conjugate to (x,y), and

k² = k² + k:

Equations (13) can be used to rewrite equation (10) in the form

(14)

In practice division by intensity is only performed if that intensity is greater than a certain threshold value (eg. 0.1% of the maximum value).

Division by k_r does not take place at the. point k_r = 0 of Fourier space; instead multiplication by zero takes place at this point, This amounts to taking the Cauchy principal value of the integral operator V_x ^~ -2

In order to quantitatively measure the phase of object it is necessary to incorporate some physical constants into the phase recovery algorithm given in Equation (14) relating to the experimental setup in use to quantify the variables k_x, k_y. This can be done by rewriting equation (14) in the following form suitable for implementation using a fast Fourier transform:

-N N where i, j e index the frequent components of F(I₊ -I_) 2 '2 where the intensity derivative d_zI(x,y) is obtained by subtracting two images T₊ and J_ separated by a distance δz , i and j are the pixel numbers on the image, and using the fact that the Fourier space step size is given by

NΔx

where the image is an N x N array of pixels of size Ax . Thus in addition to measuring the three intensity distributions it is necessary to know the pixel size Ax , defocus distance & and wavelength λ in order to make a quantitative phase measurement. All of these quantities can be readily determined: the pixel size can be determined directly for example from the CCD detector geometry (in the case of direct imaging), or by existing techniques for calibrating the transverse distance scales (in the case of an imaging system), the defocus distance can be measured directly, and the spectral distribution of the illumination can be determined either by monochromating the incident field or by analysing the spectral distribution of the radiation using existing spectroscopic methods.

An example of the phase-retrieval method implementing the solution of equation (14) can be represented by the flowchart shown in Figure 2. As shown in Figure 2 the quantitative determination of phase of a radiation wave field proceeds from a set of intensity measurements {l„} over the two spaced apart planes A and B. A measurement of central intensity I(x,y) in a selected plane parallel to and midway between the planes A and B is also obtained. The intensity measurements are performed over a defined array on each of the two planes A and B and the respective values subtracted to produce a measure of the intensity derivative. This value is multiplied by the negative of the average wave number. The data are split into two component sets and a fast Fourier transform is performed to produce the respective x and y components in the Fourier domain. A filter is then applied to the Fourier domain representations to correspond to the inversion of a differential operator reflected in the untransformed representation. An inverse Fourier transform is then performed on each of the x and y components to produce a spatial domain value from which the differential operator has been removed. A division by the central intensity I(x,y) obtained by averaging the intensity measurements over planes A and B is then performed if the intensity level is above a selected threshold level. The resultant data is again Fourier transformed and multiplied by the same filter to again correspond to the inversion of a differential operator reflected in the untransformed data. The resultant components are again inverse Fourier transformed and summed to provide a retrieved phase measurement.

It will be apparent that in general the method according to this invention can proceed from any suitable representative determination of intensity derivative or rate of change of intensity over a selected surface extending across the propagation direction and the intensity over that same surface. As will be explained in various examples these data can be obtained in a variety of ways and the method implemented to yield phase of the radiation wave field.

Rewriting equation (14) with:

Ω_x(k_x,k_v,a) = k_xk_r

Ω (k_x,k _y, ) = k_yk_r ^~

φ(x, y) = φ^w (x, y) + φ^w (^χ, y)

gives (15)

φ^ x, y) = F^ΛΩ_x (k_x,k_y, )F-^—F Ω_x(k_x,k_y , )F az

φ^ (x, y) = F Ω_y (k_x,k_y, )F—— F^ΛΩ_y (k_x,k_y ,a)F

/(χ,y) 3z

where: Φ(x, y)denotes the recovered phase,

F denotes Fourier transformation, and F^"1 denotes inverse Fourier transformation,

I(x,y) is the intensity distribution over the plane of interest, (X y) are Cartesian coordinates over the plane of interest,

(k_x,k_y) are the Fourier variables conjugate to (x,y)

k = 2π/ λ is the average wavenumber of the radiation,

λ is the average wavelength of the radiation, dl/dz is the estimate for the longitudinal intensity derivative, is the regularization parameter used to stabilize the algorithm when noise is present.

As given above, the solution to the transport of intensity equation (1) assumes a perfect imaging system. That is, there are no "aberrations" present in the optical system used to obtain the intensity data which is fed into the algorithm. Of course, no imaging system is perfect. The imperfections present in an imaging system may be quantified by a set of numbers: (16) A, ,A₂,A₃,...

which are termed aberration coefficients.

If intensity data were taken on an imperfect instrument whose imperfections were characterized by a certain set of known aberration coefficients A₁ ,A₂,A₃,... , it would be desirable if the filters Ω_x(k_xk_y, ) and Ω_y (k_xk_y, ) present in (15) could be replaced by modified filters which explicitly depend upon the aberration coefficients:

(17) Ω_x(k_xk_y, ,A ,A₂,A_i,...) and Ω_y (k_xk_y, ,A_l,A₂,A₃,...)

This would allow the imperfections of the imaging system to be eaεplicitly taken into account, leading to quantitatively correct phase retrieval using aberrated imaging systems. For the special case of a non-absorbing phase object in a radiation wave field of uniform intensity with weak (i.e. much less than 2π radians) phase variations the appropriate modified filters lead to the following functional form for the phase-retrieval algorithm:

₍1_{8 )}

where: aberratea(X/y) is the aberrated intensity measured at defocus distance δ ,

An_n are the aberration coefficients which characterize the imperfect imaging system. If a filter is defined:

(19) Ω(k_xk_y , ,A_l ,A₂,A₃,...) 1

J- 2π.δzl(k² + k²) -2∑_m ∑_n A_mnk_x'ⁿk;

Then (18) becomes:

(20) ^,y) = F-¹ΩF^-F-¹ΩF{l_aberrated(x,y)-l }

The term { I_abemued (x, y) -1 }is a measure of rate of change of intensity. I₀ intensity is a measurable constant for uniform intensity so that (20) is the same general form as (15) . Consequently the special case of aberration can be dealt with by changing the filter in the general method described above. The x and y component filters Ω^. and Ω_y are given by

(21) Ω_x = Ω_y =-^Ω

Figure 3 schematically shows an arrangement for quantitative phase amplitude microscopy. A sample is illuminated using a source of white light Kδhler illumination 15, commonly found on optical microscopes. The light is transmitted through an object 16 and collected by a microscope imaging system 17 and relayed to a CCD camera 18 or other digital imaging device having a planar imaging surface. Three images are collected: an in-focus image, I₀, and two slightly out of focus images 1₊ and I_. The defocus is obtained by suitable means such as a drive system 19 to adjust the microscope focus knob. The defocus introduced is usually quite small so that degradation in spatial resolution is minimised, although the optimal amount of defocus to use is determined by sample properties and imaging geometry such as magnification, numerical apertures, etc.

When taking the images the numerical aperture of the condenser is chosen to be less than the numerical aperture of the objective being used. If this is not the case then serious image degradation will occur, although the precise amount by which the condenser and objective numerical apertures should differ involves a tradeoff between image fidelity and spatial resolution, with the optimal difference depending on the sample properties and the optics used.

Intensity data from the collected images 1₊ and I_ are subtracted to produce a representative measure of rate of change of intensity (intensity derivative) . From this and the intensity data of collected image I₀ the method described above can be used to produce quantitative information about the magnitude of the phase shift in the image plane.

There may be cases in which it is desirable to take more than two images in order to obtain a better estimate of the intensity derivative dl/dz. A function can then be fitted to this data from which dl/dz can be computed and used in the phase determination method in place of the simple subtraction of two images normally used. It is also possible to operate this system in a reflection geometry to obtain surface topography. The principle of operation is the same, however the optics have to be folded back on themselves to form a reflection geometry - otherwise the process is identical.

For certain applications it can also be desirable to filter the light to a particular wavelength, although this is not necessary for the described imaging process as it works equally well with white light.

An implementation is shown in Figure 4. An Olympus BX-60 optical microscope 20 was equipped with a set of UMPlan metallurgical objectives and a universal condenser to provide Kδhler illumination. In order to be able to compare the results with existing imaging modes Nomarski DIG optics and a set of cover-slip corrected UplanApo objectives were also acquired for this microscope, enabling images to be taken of the same field of view using both phase retrieval and Nomarski DIC for the purposes of qualitative comparison. A 12-bit scientific grade Photometries SenSys CCD camera 21 equipped with a 1300x1035 pixel Kodak KAF-1400 CCD chip was added to the 0.5x video port on the microscope to acquire digital images of the sample.

The phase recovery technique of this embodiment of the invention requires the acquisition of defocused images. A stepper motor drive system 22 was attached to the focus knob of the microscope. This stepper motor 22 was attached to the parallel port of a 133_0fflz Pentium PC 23 also used to control the CCD camera 21, enabling full automation of the acquisition of through-focus image sequences. This data acquisition system was linked to custom software written to recover phase images from the CCD images, thereby enabling full automation of the image acquisition and data processing sequences.

The form of the differential operators used in the preferred embodiment of this invention are based on the optics of the system used to obtain the above-mentioned images. Thus, the operator takes into account the details of the optical system used to take the images. This is achieved by:

Determine the numerical aperture of the condenser, NAcondensor

Determine NA of objective, NAobjectwe and p₀bjective# the radius of the objective aperture

μ NA condenser

NA objective

(These are settings and dimensions on the microscope.)

Determine radial lateral spatial frequency, p, and longitudinal spatial frequency, η.

(These are dependent on the pixelation and position distribution of images taken in the series . )

Determine the wavelength, λ, of the radiation to be used.

The new form of the operator is

T_P +a²

where,

and

Figure 5 is a flow chart generally illustrating how T is determined by means of the above equation merely showing breakdown of the various components of the equation.

Since modifications within the spirit and scope of the invention may readily be effected by persons skilled within the art, it is to be understood that this invention is not limited to the particular embodiment described by way of example hereinabove.

Claims

1. A method of quantitative determination of a phase of a radiation wavefield including the steps of: (a) producing a representative measure of the rate of change of intensity of said radiation wave field over a selected surface extending generally across the wave field;

(c) transforming said measure of rate of change of intensity to produce a first integral transform representation and applying to said first integral transform representation a first filter corresponding to the inversion of a first differential operator reflected in said measure of rate of change of intensity to produce a first modified integral transform representation;

(d) applying an inverse of said first integral transform to said first modified integral transform representation to produce an untransformed representation;

(e) applying a correction based on said measure of intensity over said selected surface to said untransformed representation; (f) transforming the corrected untransformed representation to produce a second integral transform representation and applying to said second integral transform representation a second filter corresponding to the inversion of a second differential operator reflected in the corrected untransformed representation to produce a second modified integral transform representation;

(h) wherein at least one of the first or second differential operator has a form based on an optical system used to acquire the radiation for producing the representative measure of the rate of change of intensity of the radiation wavefield over the selected surface extending generally across the wavefield.

2. The method of claim 1 wherein both the first and second differential operators have a form based on the optical system.

3. The method of claim 2 wherein the first and second integral transforms are produced using a Fourier transform.

4. The method of claim 2 wherein the differential operators have the form:

T_P + ²

where,

and

and wherein p is radial lateral spatial frequency, the η longitudinal spatial frequency and δz is the defocus distance in the plane of the object, and

¹ N 'A condenser ξ = NA objective where NA_cmdensor and NA_objective are respectively the numerical aperture of the condensor and the objective, and p_obj is the maximum spatial frequency accepted by the objective.

5. A computer program for quantitative determination of the phase of a radiation wavefield including code to perform the method steps according to claim 1.

6. An apparatus for phase amplitude imaging of an object including: a radiation wavefield source to irradiate the object; an imaging system to focus radiation from the object to an imaging surface extending across the wavefield propagating from the object; means to produce a representative measure of radiation intensity over the imaging surface; and processing means to:

(iv) apply an inverse of said first integral transform to said first modified integral transform representation to produce an untransformed representation; (v) apply a correction based on said measure of intensity over said selected surface to said untransformed representation;

(vi) transform the corrected untransformed representation to produce a second integral transform representation;

7. The apparatus of claim 6 wherein the differential operator, in Fourier space, is

where,

and

c ε NA condenser

^~ NA objective

where NA_condensor and NA_objective are respectively the numerical aperture of the condensor and the objective, and p_obj is the maximum spatial frequency accepted by the objective. 0

8. The apparatus of claim 7 wherein both the first and second differential operators have the above form.