Image based adaptive optics through optimisation of low spatial frequencies

We present a wave front sensorless adaptive optics scheme for an incoherent imaging system. Aberration correction is performed through the optimisation of an image quality metric based upon the low spatial frequency content of the image. A sequence of images is acquired, each with a different aberration bias applied and the correction aberration is estimated from the information in this image sequence. It is shown, by representing aberrations as an expansion in Lukosz modes, that the effects of different modes can be separated. The optimisation of each mode becomes independent and can be performed as the maximisation of a quadratic function, requiring only three image measurements per mode. This efficient correction scheme is demonstrated experimentally in an incoherent transmission microscope. We show that the sensitivity to different aberration magnitudes can be tuned by changing the range of spatial frequencies used in the metric. We also explain how the optimisation scheme is related to other methods that use image sharpness metrics. © 2007 Optical Society of America OCIS codes:(010.1080) Adaptive Optics; (010.7350) Wave-front sensing; (110.4850) Optical transfer functions. References and links 1. R. K. Tyson,Principles of Adaptive Optics, Academic Press, London, 1991. 2. J. W. Hardy,Adaptive Optics for Astronomical Telescopes, Oxford University Press, 1998. 3. R. A. Gonsalves, “Phase retrieval and diversity in adaptive optics,” Opt. Engineering 21, 829-832, 1982. 4. D. R. Luke, J. V. Burke and R. G. Lyon, “Optical Wavefront Reconstruction: Theory and Numerical Methods,” SIAM Review44, 169-224, 2002. 5. R. A. Muller and A. Buffington, “Real-time correction of atmospherically degraded telescope images through image sharpening,” J. Opt. Soc. Am. 64, 1200–1210, 1974. 6. A. Buffington, F. S. Crawford, R. A. Muller, A. J. Schwemin, and R. G. Smits, “Correction of atmospheric distortion with an image-sharpening telescope,” J. Opt. Soc. Am. 67, 298–303, 1977. 7. M. A. Vorontsov, G. W. Carhart, D. V. Pruidze, J. C. Ricklin, and D. G. Voelz, “Image quality criteria for an adaptive imaging system based on statistical analysis of the speckle field,” J. Opt. Soc. Am. A 13, 1456–1466, 1996. 8. M. A. Vorontsov, G. W. Carhart, and J. C. Ricklin, “Adaptive phase-distortion correction based on parallel gradient-descent optimization,” Opt. Lett. 22, 907–909, 1997. 9. N. Doble,Image Sharpness Metrics and Search Strategies for Indirect Adaptive Optics. PhD thesis, University of Durham, United Kingdom, 2000. 10. J. R. Fienup and J. J. Miller, “Aberration correction by maximizing generalized sharpness metrics,” J. Opt. Soc. Am. A 20, 609–620, 2003. 11. L. Murray, J. C. Dainty, and E. Daly, “Wavefront correction through image sharpness maximisation,” in Proc. S.P.I.E., ‘Opto-Ireland 2005: Imaging and Vision’ 5823, 40–47, 2005. 12. M. J. Booth, “Wave front sensor-less adaptive optics: a model-based approach using sphere packings,” Opt. Express14, 1339–1352, 2006. #81604 $15.00 USD Received 29 Mar 2007; revised 11 Jun 2007; accepted 11 Jun 2007; published 14 Jun 2007 (C) 2007 OSA 25 June 2007 / Vol. 15, No. 13 / OPTICS EXPRESS 8176 13. M. J. Booth, “Wavefront sensorless adaptive optics for large aberrations,” Opt. Lett. 32, 5–7, 2007. 14. W. Lukosz, “Der Einfluß der Aberrationen auf die optische Üb rtragungsfunktion bei kleinen Orts-Frequenzen,” Optica Acta10, 1–19, 1963. 15. J. Braat, “Polynomial expansion of severely aberrated wavefronts,” J. Opt. Soc. Am. A4, 643–650, 1987. 16. V. N. Mahajan,Optical Imaging and Aberrations, Part II. Wave Diffraction Optics, SPIE, Bellingham, Wash., 2001. 17. M. J. Booth, T. Wilson, H.-B. Sun, T. Ota, and S. Kawata, “Methods for the characterisation of deformable membrane mirrors,” Appl. Opt. 44, 5131–5139, 2005. 18. W. Press, S. Teukolsky, W. Vetterling, and B. Flannery, Numerical Recipes in C, Cambridge University Press, 2nd ed., 1992. 19. J. P. Hamaker, J. D. O’Sullivan, and J. E. Noordam, “Image sharpness, Fourier optics, and redundant-spacing interferometry,” J. Opt. Soc. Am. 67, 1122–1123, 1977. 20. M. A. A. Neil, M. J. Booth, and T. Wilson, “New modal wavefront sensor: a theoretical analysis,” J. Opt. Soc. Am. A 17, 1098–1107, 2000.


Introduction
The objective of all adaptive optics systems is to reduce the wave front aberrations to an acceptable level.Normally this would involve a wave front sensor to measure the aberrations, which are in turn corrected using an adaptive element, such as a deformable mirror [1,2].In imaging systems, however, direct wave front sensing is not straightforward and wave front sensorless schemes are often employed.In certain situations, some aberration information can be extracted from a single image using phase retrieval methods; further information is obtained from two or more defocused images using the methods of phase diversity [3].These methods use iterative calculations based upon a model of the imaging process to retrieve the aberrations and the object structure.However, these calculations are not guaranteed to converge to a unique solution for arbitrary objects [4].In other wave front sensorless systems, the adaptive element is reconfigured in order to optimise a metric related to image quality.The optimisation procedure involves measurement of the metric for a number of trial correction aberrations, followed by the estimation of an improved correction aberration.This process is repeated until the image quality is considered acceptable.The number of measurements required during this process depends upon the optimisation algorithm and parameters used, the mathematical representation of the aberration, and the object structure.
Most previous work in this area has used model-free algorithms, such as simplex optimisation, conjugate gradient search or multidithering [5,6,7,8,9,10,11].These schemes have employed several optimisation metrics that are appropriate for image-based adaptive optics.However, the derivation of a constructive model based upon these metrics is complicated by the image formation process, which depends on both the aberrations and the object structure.An effective model-based adaptive optics scheme should be object independent, so the model should permit the separation of aberration and object influences on the measurements.We show that this separation is possible through the appropriate choice of optimisation metric and aberration representation.A similar approach has been demonstrated for adaptive optics in focussing systems.Using a Strehl intensity metric in conjunction with a Zernike mode aberration expansion has led to efficient schemes for the correction of small aberrations [12].Another system capable of correcting larger aberrations has been demonstrated, using a metric related to the focal spot radius and an alternative aberration expansion in Lukosz modes [13].
In this paper we describe and demonstrate an image-based adaptive optics scheme that is predominantly independent of object structure.This scheme uses the low spatial frequency content of the image as the optimisation metric but leads to correction for all spatial frequencies.The aberration is represented in terms of Lukosz modes; these modes are ideal for modelling the effects of aberrations on the imaging of low spatial frequencies [14].We describe the imaging process in terms of spectral densities and the optical transfer function.The optimisation metric g is introduced as the sum of a range of low frequencies and is related to the coefficients of the aberration expansion, {a i }.Because of this choice of aberration expansion and optimisation metric, the function g({a i }) is found to have a paraboloidal maximum that permits the use of simple maximisation algorithms.Moreover, we show that this optimisation can be performed as a sequence of independent maximisations in each aberration coefficient.The correction scheme is demonstrated for imaging in an incoherent transmission microscope.

Image formation in an incoherent imaging system
For incoherent imaging, the image I(x) is given by the convolution of the object function t(x) and the intensity point spread function (IPSF), h(x), of the system: where x is the position vector in the image plane; for clarity we have omitted magnification factors.The object is, of course, independent of any aberrations in the optical system and all aberration effects are therefore manifested in the IPSF.If instead we consider the imaging process in the frequency domain, the image Fourier transform (FT), J(m), is given by where m is the spatial frequency vector, H(m) is the optical transfer function (OTF), which is equivalent to the FT of h(x), and T (m) is the FT of t(x).In general, each of the terms in Eq. 2 is a complex quantity.In order to deal solely with real quantities, we can also describe the imaging process in terms of spectral density functions.Defining the object spectral density function as S T (m)=|T (m)| 2 and the image spectral density as S J (m)=|J(m)| 2 , then where |H(m)| is the modulation transfer function (MTF).In Eq. 3, all aberration effects are confined to the MTF.

Image spectral density for low spatial frequencies
In this section we derive approximations for S J (m) at low spatial frequencies.Expressions for heavily aberrated OTFs at low spatial frequencies have been derived using the geometrical OTF [14,15,16].We present an alternative derivation, based upon the diffraction OTF, providing equivalent expressions that are valid for all aberration magnitudes.The diffraction OTF is calculated as the autocorrelation of the effective pupil function, P(r): where r is the position vector in the pupil and P * is the complex conjugate of P.This is illustrated in Fig. 1.When the pupils are circular and aberration free, we define the pupil function P(r)=Π(r), where Π(r)=1 for |r|≤1 and zero otherwise.From this we obtain the familiar expression for the OTF as where m = |m|.A normalisation factor has been introduced so that H 0 (0)=1.The spatial frequencies are also normalised such that the cut-off of the incoherent imaging system corresponds to |m| = 2.We model phase aberrations as a function Φ(r), such that the pupil function P(r)=Π(r) exp [ jΦ(r)], where j = √ −1.The corresponding, aberrated OTF can be calculated as This is valid for all spatial frequencies.However, for small spatial frequencies, the phase term Φ(r − m) can be approximated using a Taylor series expansion and Eqn.6 can therefore be approximated by where the dot represents the scalar product, ∇ is the gradient operator and O(m 2 ) represents error terms of at least the second order in m.For small arguments, the exponential term can also be expanded as a Taylor series, so that Eq. 7 becomes The first integral is equivalent to H 0 (m).For the other integrals, the effective region of integration is defined by the overlap of the two offset pupils Π(r − m)Π(r), shown as A in Fig. 1(b).If we approximate this region by the circular pupil P, then the corresponding approximation error is at least second order in m.Equation 8 can therefore be written An approximation for the squared MTF, |H(m)| 2 , follows as where the error terms have now been omitted.The final term in Eq.We assume in the rest of this paper that these aberration modes play no role, although we note that the components are readily extracted from the phase of the image FT.If these modes are not present, Eq. 10 becomes Substitution of this expression into Eq. 3 yields an expression for the image spectral density at low spatial frequencies: In the equivalent derivation using the geometric OTF, the term H 0 (m) 2 would be replaced by unity.The corresponding error in the calculation of the MTF would be of order m, leading to significant differences compared to the diffraction OTF calculations.

Optimisation metric based upon low spatial frequencies
In this section, we derive an optimisation metric that uses the low spatial frequency content of an image.The metric permits the separation of the effects of specimen structure and aberrations.
Let g(M 1 , M 2 ) be defined as the total "energy" in all image spatial frequencies lying within the annulus for which M 1 ≤|m|≤M 2 , where M 2 is small: If we define m =(m cos ξ , m sin ξ ), then S T (m) must be a periodic function of the polar angle ξ and can be represented by its Fourier series: Note that this series has fundamental period π as the spectral density always has even symmetry about the origin, such that S T (m)=S T (−m).This permits us to calculate the first integral with respect to ξ in Eq. 13 as In order to calculate the the second integral with respect to ξ , we can first show that the pupil integral, by expanding the scalar product, becomes where χ(r) is the polar angle of the vector ∇Φ(r).The second integral is then calculated as Hence, we find that only the non-azimuthally variant component α 0 (m) and the first order terms α 1 (m) and β 1 (m) contribute to the value of g.Significant values of these first order coefficients will indicate that the object has noticeable periodicity in a predominant direction, such as a one dimensional grid.For other object structures more likely to be encountered in applications α 1 (m) and β 1 (m) are expected to be small so the corresponding terms in Eq. 17 can be neglected.Hence, we find that where and are both positive quantities if α 0 (m) = 0 in the frequency range of interest.If M 1 and M 2 are fixed and the object contains frequencies in this range, it can be seen from Eq. 18 that g is maximum if and only if |∇Φ(r)| = 0 for all r or equivalently when Φ(r) is a constant.
Although g is based only on low spatial frequencies, the optimisation process will remove all phase aberrations and hence improve imaging quality for all spatial frequencies.It is therefore appropriate to use g as an optimisation metric for aberration correction.It can also be seen from Eq. 18 that the variation of g for different aberrations can be derived entirely from the properties of the integral

Aberration expansion in Lukosz functions
In order to understand the effects of the aberration on I 1 , it is useful to represent the aberration as a combination of Lukosz functions.These functions, based upon the Zernike polynomials, were first derived by Lukosz [14] and later, independently by Braat [15].Like Zernike circle functions, the Lukosz functions are each expressed as the product of a radial polynomial and an azimuthal function and use the same dual index and numbering scheme.They can be defined as where n and m are the radial and azimuthal indices, respectively, and R m n (r) is the Zernike radial polynomial given by It is also convenient to refer to the Lukosz polynomials using a single index numbering scheme, which is explained in Appendix B. We express the aberration as a series of N where the piston, tip and tilt modes (i = 1, 2, 3 respectively) have been omitted.Using this aberration expansion, we find that each mode contributes independently to I 1 : Note that, in contrast to the derivations of Lukosz and Braat, we have chosen the normalisation of the radial polynomials B m n (r) to ensure that the weighting of each coefficient in Eq. 26 is independent of the coefficient's indices.The normalisation of B m n (r) used here is also slightly different to that employed in Reference [13].
The optimisation metric g in Eq. 18 can now be directly expressed in terms of the set of aberration coefficients, {a i },togive where for clarity we have omitted the explicit dependence on M 1 and M 2 .This approximation is accurate for small aberration amplitudes.However, for larger amplitudes it can give inaccurate, even negative values, whereas in practice g would tend to zero.A more appropriate approximation is a Lorentzian function, which is always positive, tends to zero for large aberrations, and retains an identical form to relation 27 for small aberrations: where q 2 = 1/q 0 and q 3 = q 1 /q 2 0 .This Lorentzian approximation provides a close fit to empirical measurements of g, as shown in the next Section.

Experimental investigation of the optimisation metric
The properties of the optimisation metric g({a i }) were investigated experimentally using the system shown in Fig. 2(a).The system comprised an incoherent transmission microscope with a deformable mirror (DM) and a CCD camera.For the purposes of this demonstration, the DM acted as both aberration source and correction element.A light emitting diode (LED; Lumileds, Luxeon Star/O, centre wavelength 650nm) provided incoherent illumination to a transmissive specimen which was imaged using a 150mm focal length achromatic doublet as the objective lens.An iris provided the 5mm diameter limiting aperture of the imaging system at the pupil plane of the objective.This pupil plane was imaged onto the DM (Boston Micromachines Corp., Multi-DM, 140 element, 1.5µm stroke) using a 4f system.The DM was then re-imaged through the same 4f system onto the pupil plane of the tube lens, which formed an image of the specimen on the CCD camera.
The DM was driven using a set of control signals, {c i }, where each control signal was proportional to the square root of the corresponding electrode voltage.This was found to produce linear operation over most of the DM's deflection range.In order to produce desired combinations of Lukosz modes, the control signals were obtained from Lukosz modal coefficients {a i } through a matrix multiplication: where the elements of the vectors a and c are identical to the elements of {a i } and {c i }, respectively.The matrix D provides conversion from Lukosz coefficients into Zernike coefficients (see Appendix B).The pseudo-inverse matrix B † permits the calculation of the control signals from the Zernike coefficients.This matrix was obtained using an interferometric method described in Reference [17], which also enabled flattening of the initial DM aberration figure.
In order to characterise the properties of g, we used a holographic scatterer (Physical Optics Corp.) as the transmissive specimen.An aberration-free image of the scatterer is shown in Fig. 2(b).This specimen is ideal for this characterisation as it contains all spatial frequencies within the pass band of the imaging system; this can be seen in the image spectral density (Fig. 2(c)).Figure 3(a) shows the measured variation of g with the root mean square (rms) aberration amplitude using different spatial frequency ranges.The aberration consisted of eight Lukosz modes (i = 4 to 11).The rms amplitude was calculated from the Lukosz coefficients as a = |a| (see Appendix B).Each data point shows the mean and standard deviations for an ensemble of 200 random aberrations of magnitude a.Each aberration was constructed by generating random coefficients in the range (-1,1) with uniform probability; the resulting vector was then scaled to the magnitude a.When only small spatial frequencies are used in the calculation of g, the deviation from the mean is small and the response is predominantly quadratic, as predicted by Eq. 27.When larger frequencies are also included, so that the low frequency approximations no longer hold, the value of g drops off more quickly and the deviation from the mean is more significant.However, the Lorentzian approximation is found to provide a close fit to the mean value of g for all curves.In Fig. 3(b), it can be seen that the width of the experimentally determined g response fits the theoretical prediction for low spatial frequencies.
It is important to note that the width of the g response can be tuned by changing the spatial frequency range.In other words, the use of smaller spatial frequencies in the metric permits the measurement of larger aberrations.This property presents the possibility of schemes where aberrations are corrected in a series of optimisations covering first the large magnitude aberrations (using the smallest spatial frequencies), progressing to correction of less significant aberrations using larger frequencies.

Optimisation scheme
The aberration correction process can be performed as the maximisation of g({a i }).Using relation 28, we see that this is equivalent to the minimisation of a different metric G({a i }) defined as where the approximation is valid over all aberration amplitudes.Using the Lukosz coefficients {a i } as an N-dimensional coordinate basis, it is clear from Eq. 30 that G({a i }) has a uniform paraboloidal shape in the neighbourhood of its minimum.This representation is particularly advantageous for optimisation, as the minimum of a paraboloidal function is readily found from a small number of metric evaluations.Moreover, the minimisation of G can be decomposed into a sequence of N independent one dimensional parabolic minimisations in each of the coefficients a i .In order to perform a minimisation with respect to the coefficient a k of a particular mode L k , we can express G as where can be considered a constant.As the values of q ′ 2 and q 3 are not known, the value of a k that minimises G(a k ) can be calculated from a minimum of three measurements of G.In practice, we took these three measurements by intentionally introducing different aberrations using the adaptive element.We refer to these aberrations as biases.The biases were chosen to be Φ = −bL k , Φ = 0 and Φ =+bL k , where b is a suitable bias amplitude.An image was acquired and its FT and spectral density were calculated.The appropriate range of frequency components was summed, giving the metric measurements g − , g 0 and g + respectively, and the reciprocal of each result was calculated, giving G − , G 0 and G + .The optimum correction aberration was then estimated by parabolic minimisation as [18] which is exactly equivalent to the Lorentzian maximisation of the metric g.To correct this single mode, the correction aberration Φ = a corr L k would be added to the deformable mirror.For multiple mode correction, each modal coefficient would be measured in this manner before applying the full correction aberration containing all modes.

Correction of a single mode
The correction process is illustrated in Fig. 4 for the correction of one Lukosz mode using the scatterer specimen.A suitable range of spatial frequencies and the bias amplitude were chosen based upon the curves in Fig. 3.An initial aberration was added using the DM, an image was acquired and the value of g was calculated.Positive and negative bias aberrations were added in turn and the corresponding values of g were calculated.The correction aberration was obtained using Eq.33 and the correction was applied to the DM.The final rms phase aberration was found to be 0.18, corresponding to a Strehl ratio of 0.97.

Correction of multiple modes
The correction of multiple aberration modes is illustrated in Fig. 5 for both the scatterer specimen and a US Air Force test chart.In each case, the initial aberration, a in , introduced by the DM consisted of eight Lukosz modes (i = 4 to 11) and had a total amplitude of |a in | = 11.5.Each modal coefficient was estimated in turn using a bias amplitude b = 9.8.Once all eight coefficients had been estimated, the full correction aberration, a corr , was added to the DM.We note that the unbiased measurement was identical for each modal estimation, so was only taken once.The final Strehl ratios were found to be 0.87 for the scatterer and 0.91 for the test chart.A further cycle of correction was also performed (not shown in the figure) using bias b = 4.9, giving final Strehl ratios of 0.99 for the scatterer and 0.98 for the test chart.

Accuracy of correction
We investigated the correction accuracy for various spatial frequency ranges, bias amplitudes and input aberrations using the scatterer specimen.The results are summarised in Fig. 6 and Table 1.As in the previous section, the initial aberration, a in , introduced by the DM consisted of a random combination of the eight Lukosz modes i = 4 to 11.The values of optimisation metric were obtained in a similar manner and the correction aberration, a corr , was determined.We define the aberration error to be a err = a in + a corr and the error magnitude as The approximate Strehl ratio was calculated using the Gaussian approximation [16] and related to a err using the expression derived in Appendix B:  1).The upper graphs correspond to bias b = 4.9 whereas in the lower graphs b = 1.6.(a1,a2) show the mean correction error ε; (b1,b2) show the mean Strehl ratio S. A sequence of 50 measurements was taken for each input aberration magnitude |a in | and the mean values of ε and S were calculated and plotted in Fig. 6.The data show that each of the frequency combinations can provide effective correction over a particular range, but that the range is largest when the lowest spatial frequencies are used.This property is related to the width of the g(a) curves shown in Fig. 3(a).Some of the data points in Fig. 6(a1) and (a2) lie above the ε = |a in | line, indicating that the corrected aberration is larger than the initial aberration.These points correspond to metrics using higher spatial frequencies, where we do not expect the quadratic approximations used in the derivation of Eq. 27 to hold.We note that for small input aberrations there is an offset error.However, in all cases shown here, this corresponds to a Strehl ratio of greater than 0.8, close to the diffraction limit.These results indicate that, when aberration statistics are unknown, a sensible strategy would involve choosing small spatial frequencies for an initial correction.This would be accompanied by a bias that is no larger than the half width of the response curve, as shown in Fig. 3(b).If further correction is required, this could be performed using a larger range of frequencies and a corresponding smaller bias.If the maximum expected aberration magnitude is known, then the bias could be chosen to be similar to this maximum.
The effect of additional aberration modes on the correction process was investigated by including random combinations of an extra eight modes (i = 12 to 19) in the initial aberration.The original eight modes (i = 4 to 11) were corrected in the same manner as before and ε was calculated taking into account only the modes that were corrected.The results obtained when different amounts of the additional modes were present are shown in Fig. 7.The error ε shows only a small variation as the amplitude of the additional modes is increased.This illustrates that different aberration modes can be corrected independently using this procedure.

Discussion and Conclusions
We have introduced a model-based adaptive optics scheme for correcting aberrations in an incoherent imaging system.Using an optimisation metric based upon the low spatial frequency content of the image and an aberration expansion in terms of Lukosz modes, we have been able to separate the effects of the different aberration modes.This allowed the optimisation to be performed as a sequence of independent corrections of each mode.Although only low spatial frequencies are used in the optimisation process, correction of all aberrations (aside from piston, tip and tilt) results.Consequently, imaging quality is improved for all spatial frequencies and not solely the frequencies used in the optimisation metric.The correction scheme is predominantly independent of object structure -the model is valid when the low spatial frequency components are not significantly concentrated in one orientation.This would occur, for example, if the image were dominated by a one dimensional grid-like pattern.Even if the object has this form, we expect the scheme to be robust -this has been indicated by preliminary results.Although the discussion in this paper was framed in the context of an incoherent imaging system, we expect this approach also to be valid for coherent or partially coherent systems.
The optimisation metric used in this paper can be related to image sharpness measures that and is thus equivalent to the metric g if using the spatial frequency range (M 1 , M 2 )=(0, 2).The methods described in this paper could therefore be extended for use with image sharpness metrics, obviating the need to calculate the image FT.For example, they would be directly applicable if the object spectrum were dominated by low spatial frequency components.

Fig. 1 .
Fig. 1.The calculation of the incoherent optical transfer function: (a) the circular pupil P with circumference C; (b) The geometry used in the autocorrelation calculation, showing the pupil overlap A; (c) The resulting aberration-free incoherent optical transfer function.

Fig. 2 .
Fig. 2. (a) Schematic diagram of the experimental apparatus; (b) Raw image of scatterer without aberrations; (c) Spectral density of the scatterer image (log scale) with M 1 , M 2 and incoherent cut-off frequencies marked.The horizontal and vertical lines at the edge of this image are FT artefacts arising from the sharp image boundaries.

Fig. 3 .
Fig. 3. Experimental measurement of the optimisation metric: (a) Variation of g with aberration magnitude a = |a| for the different frequency ranges given by the figures in parentheses (M 1 , M 2 ).The solid lines are Lorentzian fits to the mean values; (b) The measured and calculated half width of g(a) for different frequency ranges.The theoretical curve is valid for small frequencies.

Table 1 .
Spatial frequency ranges for the results of Fig.6.