duncansapien

Chapter 3
The Weak Phase Approximation

Introduction

Chapter one introduced the Zernike phase filtering operations which result in an image intensity proportional to the object phase. Two equations commonly appear in the explanation of the filter operation - expansion of the exponential e^if(x), equation 2.2 and the expression for the resulting spectrum 1+iF(n), equation 2.3. These equations are deemed valid when the object f(x) has a relatively small phase retardance, and together summarise what is commonly referred to as the weak phase approximation.

The first section of this chapter suggests a possible redefinition of what is meant by the `weak' phase approximation in the light of the Bessel function analysis of chapter two. If the object phase retardance f(x) is represented as a Fourier series then, it shall be shown, a useful parameter in the analysis of the spectrum is the maximum size of any Fourier coefficient a<sub>j</sub>, which in this thesis is termed the `Zernike Limit'. A refined calculation of this limit, encompassing a number of spectral and image quality parameters as useful aids, is presented in section two. Section three briefly investigates likely causes of poor image quality and obtains some unexpected results from the Bessel function program described in Chapter two.

1  Classical Explanation

The classical phase contrast filtering operations introduced by Zernike were explained at the beginning of chapter two with the assistance of a Taylor series expansion of the object so that
The Fourier Transform G(n) of the phase object was found to be
where F(n) is the Fourier Transform of f(x). This line of reasoning clearly leads to an explanation of the phase contrast imaging process, but it will now be shown that an approach from a different perspective can lead to a more complete understanding of the process.

As explained in the previous chapter, each Fourier coefficient produces a whole comb of spectral orders in the frequency plane which convolves with the combs from every other coefficient a_j, where

With regard to producing a spectrum linear with F(n), only the primary orders of each comb, having a complex amplitude of J_±1(a_j), are required. Should the size of coefficients a_j be so large as to cause appreciable ghost spectral orders (J_k(a_j) : k ³ 2) to appear then the resulting convolutions with these orders effectively destroys this linearity. Therefore, if convolutions with ghost orders are assumed to be the principle cause of spectrum (and consequently image) nonlinearity an image intensity linear in f(x) should occur after a Zernike phase contrast operation if each and every a_j is less than a certain limit - the `Zernike' limit.

Two consequences of this conjecture are immediately apparent:

1. f(x) may be as large as allowed by equation 3.3 and thus may not be small.
2. The upper limit to a_j is independent of the number of terms in the series describing f(x).

Care must be taken with this conjecture - in particular convolutions occur between the primary orders of each comb and the zero orders of every other comb (or the resultant unit comb) as illustrated in figure 3.6. Also, even though the convolutions which do occur may alter the amplitude spectrum by only a small percentage, the phase will also be altered slightly. The complex amplitude at each frequency then has a phase which is roughly `i' times F_j, plus or minus a small error term due to the convolution with other Bessel combs. The introduction of this phase error, which is likely to be random in nature, is akin to viewing the spectrum through a (small retardance) random phase plate. It is well known that image formation can withstand relatively large errors in the amplitude spectrum a_j but is highly sensitive to the phase spectrum F_j. Therefore, any calculation to determine the Zernike limit to a_j based on amplitude errors alone is likely to result in an absolute maximum figure.

1.1  Linearity and Non-Weak Phase Objects

The image intensity I(x) after phase contrast filtering at the zero of frequency may no longer obey the relation
due to the large possible size of f(x) allowed by equation 3.3 but this is not to say that the intensity distribution may no longer be linearly proportional to f(x). Recall the microscopists' view of image formation (chapter1, section 1.1.4). As the object phase retardance increases light which would have gone into the undiffracted wave to form the constant background field of the image may instead be channeled into the diffracted wave, and makes up the (larger) image field describing the form of the object f(x). This idea is illustrated in figure 3.1.

The principle of conservation of energy requires that ò_-¥^+¥I(x) dx=Constant for both images shown, and that in a loss free system¹ the constant is the same for both images.

The upper (Zernike) limit to the size of a_j may be determined from close inspection of the Bessel function graphs, or preferably by calculation. As the argument x of J_k(x) increases, so the Bessel functions increase in accordance with their order, which is to say that J_k(x) is always greater than J_k+1(x), k ³ 1, until the functions start to decrease again at relatively large x. Therefore it is only necessary to find that value x at which J₊₂(x) becomes significant in order to determine the Zernike limit. Choosing a significant value of J₊₂(x) to be 0.03, it is found J₊₂(x) ³ 0.03 when x=0.5. In order to be doubly sure of not stepping into the region where J₊₂(x) becomes significant it would be prudent to define the Zernike limit to be 0.4.²

To summarise,

1. Classically, the weak phase approximation refers to a Taylor expansion of the function describing the complex transmission of the object plane when it contains a weak phase object of phase retardance f(x). Classically, this should not exceed [(l)/10] to be called weak [7, page 490].
2. From a convolution viewpoint, it is perhaps more appropriate that the weak phase referred to is that of the Fourier coefficients describing f(x) and not f(x) itself.
3. If ghost order convolutions dominate the formation of spectral (and therefore image) non-linearities, this would imply linear imaging of the phase f(x) is possible if no Fourier coefficient of f(x) exceeds a finite limit, here termed the Zernike limit. At this stage, it is thought such a limit should not exceed the value of 0.4 radians.
4. As a consequence of (3) above, f(x) need no longer be small and may exceed the conventional value of [(l)/10].
5. The above estimate is obtained from considerations of errors in the amplitude spectrum of f(x) and is therefore likely to be an absolute maximum.

These conjectures will be tested against results obtained by computer simulation in section 3.3. Recalling the essence of the previous chapters, convolution is at heart the fundamental process at work in determining the properties of a phase object spectrum - the Taylor expansion result is completely contained within its framework but of course the opposite is not true. As such, conjectures based on investigation into convolution effects should be given more weight over the rather limited predictions of the Taylor expansion.

1.2  A Qualification

It shall now be argued that the range of Fourier coefficients a_j affects the accuracy with which the image intensity distribution approximates the object phase f(x).

This result can best be demonstrated by considering the convolution of the N=2 comb with the first unit comb. So as to reduce the amount of algebra to a minimum, an observation frequency g = 1 is taken again for which equation 3.41 of Chapter 2 shows the complex amplitude after the first convolution stage to include the terms

As mentioned in Chapter 2, the first of these terms contains the information on the spectrum of f(x) and the remaining terms are all non-linear functions of the ideal amplitude [(a₁)/2] and phase f₁ (all f_j set to zero in the above expression for compactness) at the primary order location g = +1. It is pertinent to ask the question `How do these additional cross-product terms perturb the spectrum from the ideal linear term J₊₁(a₁) J₀(a₂) i⁺¹ ? '

Considering the non-linear terms to be noise terms, the noise to signal ratio for this frequency may be defined as

for any set values of a₁ and a₂. If the phases of each Fourier coefficient are taken into account, as well as the i^m₁+m₂ terms, then the noise function is actually a phasor addition. However, increasing a_j alone serves to increase the modulus of each component phasor so the modulus of the final complex number is also increased. For the purpose of illustration only, it is deemed acceptable to define the noise to signal ratio as above. The gradient [ d/dx]J_n(a_j) is always greater than zero for n ³ 1 in the range of a_j concerned with here, and as da_j is positive the numerator of the noise to signal ratio must increase as the range of a_j increases. At the same time, the denominator must decrease as [ d/dx]J₀(a_j) is negative. Thus the noise to signal ratio must always increase if the range of coefficients increases.

The ideal spectrum of this simple phase object would only have spectral orders at primary order locations g = ±1, g = ±2 as well as the DC order. Having considered the noise to signal ratio at g = +1, what of the N=2 primary order location g = 2 ? As a_j increases not only do the ghost spectral orders increase but the linearity of J₁(a_j) with a_j decreases, and it is primarily this effect which will decrease the denominator of the noise to signal ratio at g = 3. If the noise to signal ratio of this unit comb has been increased, it would be surprising if a further convolution with this comb (to find the spectrum of a three frequency phase object) would result in a decreased noise to signal ratio of the next unit comb as a whole.³ Further, the above calculation could be performed for the g = 3 frequency position to show that increasing a₃ increases the noise to signal ratio there and so on for any number of convolutions. This argument suggests that to obtain an image of intensity linearly proportional to the object phase, not only must a_j lie below the Zernike limit but (from consideration of equation 3.6) the spectral distribution of a_j should be as flat as possible.

2  Numerical Simulation

A FORTRAN computer program was written to test the conjectures of this chapter. Two types of Fourier coefficient sets were used in the simulations, one having a_j and phase F_j selected at random and the other with a_j set to follow a Gaussian distribution with F_j again being randomly selected. The maximum value of the Fourier coefficients was incremented in steps of 0.02 and, in order to obtain statistically valid results, 50 random data sets (and therefore phase objects), were made for each set of limits to a_j.

For the first run of the program (using random a_j) twelve coefficients were chosen to make up each phase object. The minimum limit to the size of a_j was set to be 0.02, and the maximum allowed to vary in steps of 0.02 up until the Zernike limit of 0.4 was reached. The second run used a Gaussian distribution of a_j and the number of Fourier coefficients `N' selected so that 95% of the energy of the Gaussian was contained within -N < j < +N. The coefficients were defined as

and s was chosen to be 16, requiring 28 Fourier coefficients a_j to fulfil the 95% energy condition. Typical computed Gaussian amplitude spectra (where deviations from the smooth shape show up best) for three values of a_max are shown in figure 3.3.

Upon Fourier Transformation, the amplitude and phase of each spectrum was recorded and compared with that of the ideal spectrum of f(x). As the image intensity function may resemble f(x) in form but be scaled and DC offset, it is necessary upon re-transformation to find a scaling parameter which resulted in the best fitting of object phase f(x) to the shape of the image intensity. Appendix six provides a brief summary of the scaling procedure and some helpful information on the fine details of the simulation such as selecting object wavelengths to eliminate numerical phase errors in the FFT routine.

3.1  Spectrum Quality

For each spectrum in the simulation, the mean spectrum amplitude and phase deviations, D_amp and D_F, were calculated. If F_j represents the amplitude and q_j the phase of the measured spectrum (in units of p) then
where [`(a<sub>j</sub>)] is the mean of the amplitudes a<sub>j</sub> and N is the total number of coefficients. Equation 3.2 shows the object spectrum in frequency space is multiplied by `i', so that for a direct comparison between object phase and numerically recovered spectral phase data it is first necessary to subtract [(p)/2] radians from F_j. The [`(a_j)] term in the denominator of D_amp is merely a weighting, the relative errors of each frequency being weighted by the importance of that frequency compared to the mean of all amplitudes a_j, so that each term in the sum is really
The merit of such a weighting may be argued but as the simulations are primarily intended to check on image linearity, this is not a great point of contention.

The mean and standard deviation of D_amp and D_F (which are themselves mean values) are calculated for the 50 data sets found for each upper limit to a_j. In order to obtain a parameter characterising the range of a_j, the Fourier coefficients are sorted into size order and a least squares algorithm computes the best gradient m so that

where the b_I are the size ordered a_j. This is a slightly more accurate method of estimating the range of a_j than taking the parameter [(a_max-a_min)/N], especially so in the case where a_j obey a Gaussian distribution.

3.2  Ultrasonic Analogy

Strong similarities exist between the spectra of phase objects and the spectra recorded from Doppler velocimetry measurements using ultrasound. As such, simulations of ultrasonic spectra [26], [27] have been analysed in an almost identical manner to the analysis of phase object spectra in this chapter and this area is briefly reviewed here.

In Doppler velocimetry, a highly directional, very narrow bandwidth sound pulse is reflected off various surfaces (scatterers) within the medium under study and the time taken to receive the reflection is taken as proportional to the distance of the reflector from the source. If the reflector is moving then a frequency shift will occur which can be detected and used to obtain information about the velocity profile within the medium. The analogy with phase object spectra comes about due to the interaction of the reflected signal from every scatterer with the signal from every other scatterer. Due to the velocity distribution within the medium (typically blood), the sound waves reflected from each scatterer have a slightly different frequency. Further, a reflected wave may scatter again and again off many scattering centres before reaching the detector.

The signal received at the detector is the Fourier integral of the signals from these reflected waves. Suppose n scattering particles exist within the medium. The ideal Doppler spectrum would then consist of n single peaks separated from the source frequency by known amounts dn_i according to the formula

where n_s and c are the source frequency and speed of ultrasound in the medium.⁴ The reflected waves from each scatterer interfere and, due to the particle velocity distribution, result in a measured Doppler time signal [26] of
where the a_m random Rayleigh variables⁵ and F_m are uniformly distributed on the range -p,+p. M is the number of bins into which the frequency range of the spectrum has been divided - essentially the spectrum resolution. In the case of a phase object every sine wave component describing the phase retardance creates not just one but a whole spectrum of frequencies in the Fourier plane. Further, each of these frequencies interacts with the frequency combs produced by every other sinusoid of the object. In ultrasonics the word ascribed to these unwanted frequencies is speckle and its suppression is desirable in that only the underlying first order spectrum (from the initially scattered waves) contains useful information. This is analogous to the effect of the ghost orders on the image of a phase object studied in this chapter.

In particular, the simulations studied here compare the spectrum with that of the ideal spectrum which is known exactly. In Doppler ultrasonics the time average Fourier Transform of equation 3.13 is used as the goal for which all speckle reduction techniques should aim to reach, and equations identical in form to equation 3.8 of this chapter seek to quantify the variation of the spectrum from the ideal; the only difference being that in the ultrasonic case the a_j and F_j naturally vary with time whereas in this chapter they are artificially varied to simulate different phase objects.

3.3  Image Quality

Before any quantitative observations are made it is also necessary to somehow quantify the quality of the resulting image. Many possible measures of image quality exist but only two have been selected for this simulation - the cross-correlation of object phase with (scaled) image intensity, and the fidelity defect [7, pages 171-175] of the image as compared to the object. Cross-correlation is a familiar tool used to quantify differences between two functions, though the fidelity defect may not be. Briefly, if g(x) represents the image intensity of an object with amplitude transmittance f(x), the fidelity of the image is defined as
and as such is a measure of the mean difference between both functions. In the context of this chapter, f(x) is not the object transmittance but the object phase which it is hoped to produce an image of after an intermediate spatial filtering process. In order to gain a single figure for each image quality parameter, only the central ordinate of the cross-correlation was calculated, which reduces to
Figure 3.3 shows both the mean amplitude and phase deviations together with the cross-correlation and fidelity defect results obtained from the first run of the program, and figure 3.4 for the second run. Graphs measuring the same parameter are drawn to the same scale for both runs for the program so that visual comparison may readily be made.

3.4  Observations

Firstly it is clear that both spectral integrity and image linearity deteriorate as the mean amplitude of Fourier coefficient increases. This implies one of two things - either the ghost orders are still making a significant contribution to the spectrum or convolutions with the primary orders cause more harm than previously imagined. In either case one must conclude that detrimental convolution effects are still occuring. The conjecture that the spectrum deviation increases as the range of Fourier coefficients increases is verified in this simulation even at extremely small values of a_max.

Tables 3.1 and 3.2 contain the results of three regions of interest from the first and second run of the program respectively. Examining the image quality results obtained from the first run (random a_j) observe that the mean cross-correlation between object phase and image intensity falls to @ 88±8% when a_max=0.4 (corresponding to [`(a_j)] @ 0.2). It is the opinion of the author that this value lies in the lower range of what might be considered an acceptable figure for this parameter. Reducing the value of a_max to 0.2, it is observed, results in a much improved cross-correlation and fidelity defect which are better still when a_max falls to 0.1.

For a given a_max, the second run shows that objects with spectra obeying a Gaussian distribution of amplitude generally suffer more severe spectrum and image quality degradation than do objects with random spectral amplitude profiles. For the random spectral profile the higher frequencies will, on average, contain as much energy as the lower frequencies and one might expect the detrimental convolution effects to cause roughly equal amounts of spectral deviance from the ideal over the whole frequency range. However, by its very nature the Gaussian amplitude distribution has less energy in the higher frequencies where convolution effects might cause relatively more harm. As it is the high spatial frequency components which are generally responsible for the image acuity, one might argue that the Gaussian spectral profile should be more sensitive to the convolution process than the random spectral profile. For the Gaussian profile, figure 3.4 shows that the increase of image fidelity defect is almost exponential in nature, and the fall off of cross-correlation as a_max increases is rapid indeed. As for the first run however, an upper limit of a_max @ 0.1 is observed to produce excellent cross-correlation and fidelity defect results.