- Email: [email protected]
Constrained least-squares estimation in deconvolution from wave-front sensing
15 May 1998
Optics Communications 151 Ž1998. 93–100
Full length article
Constrained least-squares estimation in deconvolution from wave-front sensing S.D. Ford a
a,1
, B.M. Welsh
a,2
, M.C. Roggemann
b
Department of Electrical and Computer Engineering, Air Force Institute of Technology, 2950 P Street, Bldg 640, Wright-Patterson Air Force Base, OH 45433-7765, USA b Department of Electrical Engineering, Michigan Technological UniÕersity, 1400 Townsend DriÕe, 121 EERC Bldg, Houghton, MI 49931-1295, USA Received 5 November 1997; revised 2 February 1998; accepted 4 February 1998
Abstract We address the optimal processing of astronomical images using the deconvolution from wave-front sensing technique ŽDWFS.. A constrained least-squares ŽCLS. solution which incorporates ensemble average DWFS data is derived using Lagrange minimization. The new estimator requires DWFS data, noise statistics, OTF statistics, and a constraint. The constraint can be chosen such that the algorithm selects a conventional regularization constant automatically. No ad hoc parameter tuning is necessary. The algorithm uses an iterative Newton-Raphson minimization to determine the optimal Lagrange multiplier. Computer simulation of a 1 m telescope imaging through atmospheric turbulence is used to test the estimation scheme. CLS object estimates are compared with those processed via manual tuning of the regularization constant. The CLS algorithm provides images with comparable resolution and is computationally inexpensive, converging to a solution in less than 10 iterations. q 1998 Elsevier Science B.V. All rights reserved. Keywords: Atmospheric turbulence; Constrained least-squares; Deconvolution from wavefront sensing
1. Introduction Atmospheric turbulence and its limits on the resolution of ground-based telescopes have been well documented w1–3x. Turbulent atmospheric motion causes spatial and temporal fluctuations in the index of refraction of air. These index of refraction fluctuations impose random phase aberrations on the incoming light which results in a general broadening of the atmospheric-optical system point
1 Present address: Air Force Research Laboratory ŽAFRLrDEBS., Bldg 422, 3550 Aberdeen Avenue SE, Kirtland AFB, NM 87117-5776, USA. 2 Present address: Mission Research Corporation, 3975 Research Boulevard, Dayton, OH 45431-2108, USA. E-mail: [email protected]
spread function ŽPSF. w1x. Three general classes of techniques are available to overcome turbulence effects w3x: Ži. pre-detection processing also known as adaptive optics ŽAO. w4x which senses and corrects for turbulence-induced aberrations before the light is collected; Žii. post-detection processing, such as speckle imaging w5–7x, which corrects for turbulence after the light is collected by the detector; and Žiii. hybrid methods w8–10x which combine both postdetection processing and AO hardware. One notable hybrid technique is deconvolution from wave-front sensing ŽDWFS. w10–18x. In DWFS, a wave-front sensor ŽWFS. is used to measure the pupil plane phase aberration associated with each short exposure image. This phase measurement can be used to estimate the Fourier transform of the PSF or optical transfer function ŽOTF.. The short exposure image and OTF estimate can then be used to estimate the true object irradiance via a deconvolution filter. In this
0030-4018r98r$19.00 q 1998 Elsevier Science B.V. All rights reserved. PII S 0 0 3 0 - 4 0 1 8 Ž 9 8 . 0 0 0 7 4 - 1
94
S.D. Ford et al.r Optics Communications 151 (1998) 93–100
paper, we present an optimal technique for processing noisy DWFS data based on constrained least-squares ŽCLS. estimation. The DWFS technique was first proposed by Fontanella w10x, extended by Fried w12x, and later developed by Primot et al. w11x. Primot et al. also conducted the first performance analysis of DWFS w11x. A variant of their estimator was later validated on astronomical data w13x. Welsh and Von Niederhausern further investigated DWFS performance by incorporating an optimal wave-front phase estimate w14x. Roggemann et al. w15x showed that the Primot estimator was biased and suggested a related unbiased estimator. Roggemann and Welsh also derived a signal-tonoise ratio ŽSNR. expression for DWFS w17x and conducted further comparison with speckle imaging and traditional linear deconvolution w18x. This body of research has greatly extended our understanding of DWFS performance. However, the traditional estimator given by Primot et al. w11x provides unacceptable noise amplification when processing photon-limited image ensembles. The standard solution is to add a parametric regularization constant w13x or SNR term w15x to the denominator of the traditional estimator. These approaches are analogous to CLS w20x and the parametric Wiener filter w21x, respectively. The regularization is adjusted by the user based on the perceived image quality. Thus, the resultant object estimate does not satisfy any mathematical optimality criterion. We propose a modified CLS estimation scheme which provides optimal processing of noisy DWFS data. The algorithm selects the value of the regularization constant which is consistent with the ensemble average data and a constraint. This problem is solved using the Lagrange multiplier technique w20,21x. A closed form solution for the object estimate is obtained which is analogous to the traditional DWFS estimator w13,15x with the Lagrange multiplier serving as a regularization constant. An iterative approach based on a Newton-Raphson minimization w19,22x is used to find the optimal Lagrange multiplier. The iteration incorporates the statistics of both the OTF and noise. The results of this study show that CLS estimation provides high quality processing of noisy DWFS data automatically. No ad hoc tuning of the regularization constant is necessary. CLS object estimates are compared with those processed manually. In all cases, the new technique provides images with comparable resolution. In addition, the algorithm is computationally inexpensive, converging to a solution in less than 10 iterations. The remainder of this paper is organized as follows. An image degradation model is presented in Section 2 to include the effects of shot noise and detector read noise. Section 3 reviews the traditional DWFS estimator. In Section 4, we present a CLS algorithm which incorporates a new constraint based on the ensemble average DWFS data. Section 5 presents sample results and Section 6 provides a brief conclusion.
2. Image model The standard model for a linear shift-invariant imaging system can be written as w23x i Ž x , y . s hŽ x , y .) o Ž x , y . ,
Ž1.
where iŽ x, y . is the image, hŽ x, y . is the impulse response or PSF, oŽ x, y . is the object irradiance, and ) denotes two-dimensional convolution. Taking the Fourier transform of both sides of Eq. Ž1. yields I Ž u,Õ . s H Ž u,Õ . O Ž u,Õ . ,
Ž2.
where a capital letter denotes the Fourier transform of an associated lower case quantity and H Ž u,Õ . is the Fourier transform of hŽ x, y ., also known as the OTF. In this paper, the OTF is a random quantity representing the combined effect of atmospheric turbulence and fixed telescope aberrations. Atmospheric turbulence causes a random phase aberration in the telescope pupil. The corresponding OTF is given by a normalized autocorrelation of the generalized pupil function and is a random process due to the presence of the turbulence-induced phase aberrations w23x. The OTF for a diffraction-limited full aperture is a deterministic tapered low-pass filter. Random aberrations in the telescope pupil further degrade performance resulting in a more attenuated OTF, especially at high spatial frequencies w23x. The result is a significant reduction in resolution when compared to the diffraction-limit. For large uncompensated telescopes collecting long exposure images, the resolution is approximately equal to a telescope of diameter r 0 , where r 0 is the Fried parameter which reflects turbulence strength w2x. Since typical values of r 0 range from less than 10 cm to 20–30 cm, resolution is usually limited by atmospheric turbulence not the optical system. To further complicate matters, the image iŽ x, y . is also degraded during the detection process. CCD detectors, which are often used to collect DWFS data, exhibit noise effects due to the random arrival time and location of photoevents. This noise is known as shot or photon noise and is governed by Poisson statistics w2x. CCD detectors also suffer from read noise effects denoted by the quantity n r Ž x, y .. In this paper, we assume the read noise is zero mean, Gaussian distributed, and spatially uncorrelated with uniform variance sr2. Given these assumptions, a detected image dŽ x, y . can be modeled using w3x K
dŽ x , y. s
Ý d Ž x y x k , y y yk . q n r Ž x , y . ,
Ž3.
ks 1
where K is the number of photoevents in the image, Ž x k , y k . is the location of the k th photoevent, and d denotes the Dirac delta function. In this model, K, Ž x k , y k ., and n r Ž x, y . are random variables. Since detector read noise is a form of signal independent noise, n r Ž x, y . is independent of K, Ž x k , y k ., and H Ž u,Õ ..
S.D. Ford et al.r Optics Communications 151 (1998) 93–100
To investigate the modified CLS formulation, a vectormatrix expression is needed for Eq. Ž3.. The new expression includes an additive noise vector n which incorporates both shot and read noise effects based on writing the total noise as nŽ x, y . s dŽ x, y . y iŽ x, y .. Thus, Eq. Ž3. can now be written as d s Ho q n,
Ž4.
where d and o are P length vector versions of the functions dŽ x, y . and oŽ x, y ., respectively. P is the number of pixels in the detector array. The matrix H is a P = P block-circulant matrix representing the shift-invariant PSF. H and o are properly ordered to perform the discrete convolution of Eq. Ž1. w20,21x.
95
4. Modified constrained least-squares formulation Our CLS estimation problem can be stated as follows: ‘‘Given the DWFS-derived estimate ² H˜i) Ž u,Õ . Di Ž u,Õ .:, the ensemble average magnitude-squared OTF estimate ²< H˜i Ž u,Õ .< 2 :, and statistical models for the OTF and noise, find the CLS optimal estimate of the object that caused the data.’’ To accomplish this task, consider finding an object estimate oˆ that minimizes ²< < H˜i Coˆ < < 2 : ,
Ž7.
subject to the constraint < <² H˜iTd i : y ² H˜iTH˜i : oˆ < < 2 s E < < H T n < < 2 4 ,
3. Traditional estimator Raw DWFS data consists of an ensemble of short exposure images and corresponding estimates of the wave-front phase from a WFS. The phase estimates f˜ i Ž x, y . can be used to generate an estimate of the OTF H˜i Ž u,Õ . via a normalized autocorrelation of the generalized pupil function where the subscript i refers to the ith realization. The DWFS estimator suggested by Primot et al. w11x can be written as Oˆ Ž u,Õ . s
² H˜i) Ž u,Õ . Ii Ž u,Õ . : 2
¦ H˜ Ž u,Õ . ;
,
Ž5.
i
where ² : denotes an ensemble average. Eq. Ž5. does not provide acceptable reconstructions due to residual noise in the OTF estimation process. When detector noise is present, Ii Ž u,Õ . is not available and Di Ž u,Õ . must be substituted in Eq. Ž5.. The standard solution is to incorporate a manually selected regularization constant e in the denominator of Eq. Ž5. which gives a new estimator w13,14x Oˆ Ž u,Õ . s
² H˜i) Ž u,Õ . Di Ž u,Õ . : 2
¦ H˜ Ž u,Õ . ;q e
.
Ž6.
i
The regularization constant e is adjusted manually by the user based on the perceived image quality. In the next section, we investigate an algorithm which replaces manual selection of e with an automatic algorithm based on CLS estimation. The algorithm takes advantage of noise reduction through ensemble averaging by directly incorporating ² H˜i) Ž u , Õ . D i Ž u , Õ .: and ²< H˜i Ž u,Õ .< 2 :. Thus, the objective function used to derive this estimator takes on an unfamiliar form when compared to more traditional CLS applications in image processing w19–21x.
Ž8.
where C is a block-circulant constraint matrix which enforces prior knowledge of the true object such as smoothness or support, H˜i is the block-circulant estimated PSF matrix for the ith realization, oˆ is an ordered vector containing the object estimate, and < < < < 2 is the normsquared operator. Notice that this problem statement is different from the form seen in previous CLS image p ro cessin g ap p licatio n s w 1 9 – 2 1 x . H ere, th e ² H˜i) Ž u,Õ . Di Ž u,Õ .: and ²< H˜i Ž u,Õ .< 2 : quantities are incorporated in the constraint function directly. Also, the mean-square error associated with the measurement is now E< < H T n < < 2 4 instead of the norm-squared of the noise E< < n < < 2 4 used previously w19–21x. Even though Eqs. Ž7. and Ž8. are unfamiliar, we will show that the associated Fourier domain filter has the same general form as Eq. Ž6.. The appropriate objective function J Ž oˆ . for the Lagrange minimization can be written as J Ž oˆ . s oˆT C T² H˜iTH˜i : Coˆ q l < <² H˜iTd i : y² H˜iTH˜i : oˆ < < 2 y E < < H T n < < 2 4 ,
4
Ž9.
where l is a Lagrange multiplier w20,21x. We can minimize the objective function in Eq. Ž9. by differentiating J Ž oˆ . with respect to o, ˆ setting the derivative equal to zero, and solving for o. ˆ Applying the appropriate vector-matrix identities w24x to take the derivatives in Eq. Ž9. yields y1
oˆ s ² H˜iT H˜i : ² H˜iT H˜i : q g C T² H˜iT H˜i : C
ž
=² H˜iTH˜i : ² H˜iTd i : ,
/ Ž 10.
where g s 1rl. It is not convenient to evaluate Eq. Ž10. directly since C and H˜i are large matrices for standardsized detector arrays. However, these matrices are block circulant. Block circulant matrices are diagonalized by the discrete Fourier transform allowing the transform domain
S.D. Ford et al.r Optics Communications 151 (1998) 93–100
96
Rewriting the E< < H T n < < 2 4 term in the Fourier domain using Parseval’s theorem w21x yields E < < H T n < < 2 4 sE
½Ý
H ) Ž u,Õ . N Ž u,Õ .
u ,Õ
sE
½Ý ½Ý
2
5
H ) Ž u,Õ . Ž I Ž u,Õ . y D Ž u,Õ . .
2
u ,Õ
sE
5
H ) Ž u,Õ . Ž H Ž u,Õ . O Ž u,Õ . y D Ž u,Õ . .
2
u ,Õ
5
.
Ž 13. Substituting the Fourier transform of the noise model given in Eq. Ž3. into Eq. Ž13. introduces the noise statistics such that E < < H T n < < 2 4 sE Fig. 1. Satellite computer rendering used to test CLS algorithm performance. Negative image shown for clarity.
½
Ý
ž
H ) Ž u,Õ . KH Ž u,Õ . On Ž u,Õ .
u ,Õ
2
K
y
Ý exp wyj2p Ž ux k q Õyk . x y Nr Ž u,Õ . ks 1
equivalent of Eq. Ž10. to provide a simpler scalar form w20,21x. Cancelling ²< H˜i Ž u,Õ .< 2 : terms yields Oˆ Ž u,Õ . s
² H˜i) Ž u,Õ . Di Ž u,Õ . : 2
¦ H˜ Ž u,Õ . ;q g C Ž u,Õ .
2
,
Ž 11.
i
where C Ž u,Õ . is the Fourier transform of the constraint cŽ x, y . associated with the matrix C. For all results presented in Section 5, C is the identity matrix and C Ž u,Õ . s 1 for all spatial frequencies Ž u,Õ .. Thus, Eq. Ž11. is identical to Eq. Ž6.. In this case, our iterative CLS estimation scheme provides automatic selection of the regularization constant. However, Eq. Ž11. will accommodate a general Fourier domain constraint function which is tailored to prior knowledge or a specific application. Instead of tuning g manually, we use a Newton-Raphson technique to find the parameter value which minimizes the objective function J Ž oˆ . given in Eq. Ž9.. The most efficient iteration scheme is on the Lagrange multiplier l not its inverse g w19x. To implement the Newton-Raphson technique, we require the derivative of J Ž oˆ . with respect to l. This derivative is straightforward and given by
E J Ž oˆ . El
s < <² H˜iTd i : y ² H˜iTH˜i : oˆ < < 2 y E < < H T n < < 2 4 .
Ž 12.
The first norm-squared term in Eq. Ž12. is a function of the ensemble average data and the current object estimate. It can be computed conveniently using the Fourier domain quantities ² H˜i) Ž u,Õ . Di Ž u,Õ .: and ²< H˜i Ž u,Õ .< 2 :. However, the second norm-squared term E< < H T n < < 2 4 is not a function of DWFS data and must be derived. Here, we will make use of the noise model given in Section 2.
/5
,
Ž 14.
where the normalized object spectrum OnŽ u,Õ . s O Ž u,Õ .r K has been introduced to the above expression and K is the average number of photoevents per image. Now the linearity property allows us to move the expectation inside the summation over spatial frequencies Ž u,Õ . and evaluate one term of Eq. Ž14.. Standard techniques for evaluating expectations of doubly stochastic Poisson random processes use nested conditional expectations over the random quantities Ž x k , y k ., K, and H Ž u,Õ . w3x. After cancelation of terms, the result for a single spatial frequency is E < H ) Ž u,Õ . N Ž u,Õ . < 2 4 s Ž Kq Psr2 . E < H Ž u,Õ . < 2 4 .
Ž 15. Summing the result in Eq. Ž15. over all spatial frequencies gives the final result E < < H T n < < 2 4 s Ž Kq Psr2 . Ý E < H Ž u,Õ . < 2 4 .
Ž 16 .
u ,Õ
If the OTF statistics are unavailable in Eq. Ž16., the OTF estimate data can be substituted which yields E < < H T n < < 2 4 f Ž Kq Psr2 . Ý ² < H˜i Ž u,Õ . < 2 :.
Ž 17 .
u ,Õ
The result given in Eq. Ž17. is used in Eq. Ž12. to find the current derivative of J Ž oˆ . with respect to the Lagrange multiplier l. This derivative is then used in the NewtonTable 1 Average photoevents per integration time mn
KW
K
4 6 8
3234 518 83
166469 26635 4262
S.D. Ford et al.r Optics Communications 151 (1998) 93–100
97
Fig. 2. CLS algorithm comparison with manual parameter selection, r 0 s 20 cm, mn s q4, sr s 15 electrons per pixel. Ža. Single short exposure image. Žb. CLS algorithm estimate, g s 0.0015, 4 iterations. Žc., Žd., Že., Žf. Manual parameter selection, e s 0.1, 0.01, 0.0001, and 0.00001, respectively.
98
S.D. Ford et al.r Optics Communications 151 (1998) 93–100
Fig. 3. CLS algorithm output versus object brightness, r 0 s 10 cm, sr s 15 electrons per pixel. Ža. Detected data, mn s q4. Žb. CLS estimate, mn s q4, g s 0.00027, 4 iterations. Žc. Detected data, mn s q6. Žd. CLS estimate, mn s q6, g s 0.0093, 5 iterations. Že. Detected data, mn s q8. Žd. CLS estimate, mn s q8, g s 0.8180, 7 iterations.
S.D. Ford et al.r Optics Communications 151 (1998) 93–100
Raphson iteration to update the value of l and generate an object estimate using Eq. Ž10. or its Fourier domain equivalent Eq. Ž11.. The iteration continues until the object estimate meets a pre-determined stopping criterion. A variety of criteria can be used to stop the Newton-Raphson iteration w22x. Since our goal is to minimize the objective function given in Eq. Ž9., we choose to stop the iteration when J Ž oˆ . is sufficiently small. All CLS images were processed until J Ž oˆ . F 0.001. In the next section, we show the results of using this technique on simulated DWFS data.
5. Sample results The sample results shown in this section are based on a computer rendering of a typical satellite as shown in Fig. 1. The detector array was 256 = 256 elements for a total of P s 65536 pixels. The satellite was 10 m in length and in low earth orbit at a range of 500 km. The diameter of the telescope was 1 m with both imaging and wavefront sensor wavelengths set at ls 500 nm. The telescope fractional bandwidth was "5%. The object was assumed to have the same spectral signature as the sun. The WFS subaperture size was 10 cm for a total of 60 subapertures within the entrance pupil. Atmospheric turbulence and detector noise effects were modeled using an existing computer simulation w3x. For each data realization, the simulation created a random phase screen with the appropriate turbulence statistics w25x, calculated the true OTF, and formed a detected image realization d i Ž x, y .. At the same time, a WFS model was used to generate a phase estimate f˜ i Ž x, y . which was then used to compute the estimated OTF. Finally, the required quantities H˜i) Ž u,Õ . Di Ž u,Õ . and < H˜i Ž u,Õ .< 2 were accumulated and the process repeated 150 times to generate the ensemble averages. The simulation incorporates an intensity splitter which sends 40% of the photons to the image plane and 60% to the WFS. Integration times of 10 ms were used for both WFS and imaging system. We modeled the object brightness associated with visual magnitudes mn s q4, q6, and q8. The resulting photoevents per integration time for these cases are presented in Table 1, where KW is the average number of photoevents across an individual WFS subaperture. Fig. 2 displays the detected short exposure image and CLS algorithm outputs associated with excellent seeing conditions Ž r 0 s 20 cm., a relatively bright satellite object Ž mn s q4., and detector read noise representing a high quality CCD detector Ž sr s 15 electrons per pixel.. A single short exposure image is given in Ža. to illustrate shot noise and detector read noise effects. Clearly, this image is degraded by the detection process. After optimal CLS processing of the DWFS data as shown in Žb., image resolution is greatly enhanced. In this case, the CLS
99
algorithm automatically selected the regularization parameter g s 0.0015 in 4 iterations. In Žc., Žd., Že., and Žf. the regularization parameter e associated with the traditional DWFS estimator given in Eq. Ž6. was selected manually. For images Žc. and Žd., e ) g . For images Že. and Žf., e - g . Notice that the quality of the CLS image in Žb. is as good or better than the manual images. Thus, Fig. 2 illustrates that the CLS algorithm selected a reasonable value for g in this case. Not only does shot noise degrade the short exposure data, it also restricts WFS accuracy. Without a sufficiently accurate OTF estimate, DWFS performance is severely degraded. To illustrate this limitation while applying the CLS algorithm, consider Fig. 3. Here, r 0 s 10 cm and detector read noise remains unchanged from the previous cases. Images Ža., Žc., and Že. show single short exposure realizations for object brightness cases mn s q4,q 6, and q8, respectively. In Žc. and Že., light level is so low, detector read noise is clearly dominant upon visual inspection of the images. Images Žb., Žd., and Žf. give the corresponding CLS estimates where g s 0.00027, 0.0093, and 0.8180, respectively. As object brightness decreases, the output images are more blurred. This effect is consistent with deconvolution using a poor quality estimate of the OTF and noise effects driving the CLS regularization parameter selection to a value near unity.
6. Conclusion We investigated a CLS estimator that incorporates noisy DWFS data, noise statistics, and OTF statistics. This estimator selects the regularization parameter automatically. No ad hoc tuning is necessary to reduce high spatial frequency noise effects in the DWFS image. The CLS estimator uses a Newton-Raphson iteration to select a Lagrange multiplier which minimizes an objective function. The objective function uses ensemble average data directly which aids in noise suppression. Our simulation results show that the new algorithm produces DWFS images comparable in quality to manual regularization with minimal computational expense. Finally, the CLS estimator was derived to incorporate a general Fourier domain constraint. Thus, other constraint functions can be used based on the specifics of the application.
Acknowledgements This work was supported by the Air Force Office of Scientific Research, Bolling AFB, Maryland and the Air Force Maui Optical StationrMaui High Performance Computer Center Research and Development Consortium in Kihei, Maui, Hawaii.
100
S.D. Ford et al.r Optics Communications 151 (1998) 93–100
References w1x F. Roddier, The Effects of Atmospheric Turbulence in Optical Astronomy, in: E. Wolf ŽEd.., Progress in Optics-XIX, North-Holland, New York, 1981. w2x J.W. Goodman, Statistical Optics, Wiley, New York, 1985. w3x M.C. Roggemann, B.M. Welsh, Imaging Through Turbulence, CRC Press, Boca Raton, Florida, 1996. w4x J.H. Hardy, Proc. IEEE 66 Ž1978. 651. w5x A. Labeyrie, Astron. Astrophys. 6 Ž1970. 85. w6x K.T. Knox, B.J. Thompson, Astrophys. J. 193 Ž1974. 45. w7x A.W. Lohmann, G. Weigelt, B. Wirnitzer, Appl. Optics 22 Ž1983. 4028. w8x M.C. Roggemann, D.W. Tyler, M.F. Bilmont, Appl. Optics 31 Ž1992. 7429. w9x M.C. Roggemann, C.A. Stoudt, B.M. Welsh, Opt. Eng. 33 Ž1994. 3254. w10x J. Fontanella, J. Opt. ŽFrance. 16 Ž1985. 257. w11x J. Primot, G. Rousset, J.C. Fontanella, J. Opt. Soc. Am. A 7 Ž1990. 1598. w12x D.L. Fried, in: Proc. 828 SPIE Conference on Digital Image Recovery and Synthesis, 1987, pp. 127–133. w13x J.D. Gonglewski, D.G. Voelz, J.S. Fender, D.C. Dayton, B.K. Spielbusch, R.E. Pierson, Appl. Optics 29 Ž1990. 4527.
w14x B.M. Welsh, R.N. VonNiederhausern, Appl. Optics 32 Ž1993. 5071. w15x M.C. Roggemann, B.M. Welsh, J. Devey, Appl. Optics 33 Ž1994. 5754. w16x T.J. Schulz, in: Proc. 2029 SPIE Conference on Digital Image Recovery and Synthesis II, 1993, pp. 3111–320. w17x M.C. Roggemann, B.M. Welsh, Appl. Optics 33 Ž1994. 5400. w18x B.M. Welsh, M.C. Roggemann, Appl. Optics 34 Ž1995. 2111. w19x M.C. Roggemann, D.W. Tyler, Appl. Optics 36 Ž1997. 2360. w20x B.R. Hunt, IEEE Trans. Comput. C 22 Ž1973. 805. w21x R.C. Gonzales, R.E. Woods, Digital Image Processing, Addison-Wesley, Reading, Massachusetts, 1992. w22x E. Isaacson, H.B. Keller, Analysis of Numerical Methods, Dover Publications, New York, 1994. w23x J.W. Goodman, Introduction to Fourier Optics, McGraw-Hill, New York, 1968. w24x J.L. Melsa, D.L. Cohn, Decision and Estimation Theory, McGraw-Hill, New York, 1978. w25x B.M. Welsh, A Fourier Series based atmospheric phase screen generator for simulating nonisoplanatic geometries and temporal evolution, submitted to Appl. Optics, February 1997, unpublished.
Optics Communications 151 Ž1998. 93–100
Full length article
Constrained least-squares estimation in deconvolution from wave-front sensing S.D. Ford a
a,1
, B.M. Welsh
a,2
, M.C. Roggemann
b
Department of Electrical and Computer Engineering, Air Force Institute of Technology, 2950 P Street, Bldg 640, Wright-Patterson Air Force Base, OH 45433-7765, USA b Department of Electrical Engineering, Michigan Technological UniÕersity, 1400 Townsend DriÕe, 121 EERC Bldg, Houghton, MI 49931-1295, USA Received 5 November 1997; revised 2 February 1998; accepted 4 February 1998
Abstract We address the optimal processing of astronomical images using the deconvolution from wave-front sensing technique ŽDWFS.. A constrained least-squares ŽCLS. solution which incorporates ensemble average DWFS data is derived using Lagrange minimization. The new estimator requires DWFS data, noise statistics, OTF statistics, and a constraint. The constraint can be chosen such that the algorithm selects a conventional regularization constant automatically. No ad hoc parameter tuning is necessary. The algorithm uses an iterative Newton-Raphson minimization to determine the optimal Lagrange multiplier. Computer simulation of a 1 m telescope imaging through atmospheric turbulence is used to test the estimation scheme. CLS object estimates are compared with those processed via manual tuning of the regularization constant. The CLS algorithm provides images with comparable resolution and is computationally inexpensive, converging to a solution in less than 10 iterations. q 1998 Elsevier Science B.V. All rights reserved. Keywords: Atmospheric turbulence; Constrained least-squares; Deconvolution from wavefront sensing
1. Introduction Atmospheric turbulence and its limits on the resolution of ground-based telescopes have been well documented w1–3x. Turbulent atmospheric motion causes spatial and temporal fluctuations in the index of refraction of air. These index of refraction fluctuations impose random phase aberrations on the incoming light which results in a general broadening of the atmospheric-optical system point
1 Present address: Air Force Research Laboratory ŽAFRLrDEBS., Bldg 422, 3550 Aberdeen Avenue SE, Kirtland AFB, NM 87117-5776, USA. 2 Present address: Mission Research Corporation, 3975 Research Boulevard, Dayton, OH 45431-2108, USA. E-mail: [email protected]
spread function ŽPSF. w1x. Three general classes of techniques are available to overcome turbulence effects w3x: Ži. pre-detection processing also known as adaptive optics ŽAO. w4x which senses and corrects for turbulence-induced aberrations before the light is collected; Žii. post-detection processing, such as speckle imaging w5–7x, which corrects for turbulence after the light is collected by the detector; and Žiii. hybrid methods w8–10x which combine both postdetection processing and AO hardware. One notable hybrid technique is deconvolution from wave-front sensing ŽDWFS. w10–18x. In DWFS, a wave-front sensor ŽWFS. is used to measure the pupil plane phase aberration associated with each short exposure image. This phase measurement can be used to estimate the Fourier transform of the PSF or optical transfer function ŽOTF.. The short exposure image and OTF estimate can then be used to estimate the true object irradiance via a deconvolution filter. In this
0030-4018r98r$19.00 q 1998 Elsevier Science B.V. All rights reserved. PII S 0 0 3 0 - 4 0 1 8 Ž 9 8 . 0 0 0 7 4 - 1
94
S.D. Ford et al.r Optics Communications 151 (1998) 93–100
paper, we present an optimal technique for processing noisy DWFS data based on constrained least-squares ŽCLS. estimation. The DWFS technique was first proposed by Fontanella w10x, extended by Fried w12x, and later developed by Primot et al. w11x. Primot et al. also conducted the first performance analysis of DWFS w11x. A variant of their estimator was later validated on astronomical data w13x. Welsh and Von Niederhausern further investigated DWFS performance by incorporating an optimal wave-front phase estimate w14x. Roggemann et al. w15x showed that the Primot estimator was biased and suggested a related unbiased estimator. Roggemann and Welsh also derived a signal-tonoise ratio ŽSNR. expression for DWFS w17x and conducted further comparison with speckle imaging and traditional linear deconvolution w18x. This body of research has greatly extended our understanding of DWFS performance. However, the traditional estimator given by Primot et al. w11x provides unacceptable noise amplification when processing photon-limited image ensembles. The standard solution is to add a parametric regularization constant w13x or SNR term w15x to the denominator of the traditional estimator. These approaches are analogous to CLS w20x and the parametric Wiener filter w21x, respectively. The regularization is adjusted by the user based on the perceived image quality. Thus, the resultant object estimate does not satisfy any mathematical optimality criterion. We propose a modified CLS estimation scheme which provides optimal processing of noisy DWFS data. The algorithm selects the value of the regularization constant which is consistent with the ensemble average data and a constraint. This problem is solved using the Lagrange multiplier technique w20,21x. A closed form solution for the object estimate is obtained which is analogous to the traditional DWFS estimator w13,15x with the Lagrange multiplier serving as a regularization constant. An iterative approach based on a Newton-Raphson minimization w19,22x is used to find the optimal Lagrange multiplier. The iteration incorporates the statistics of both the OTF and noise. The results of this study show that CLS estimation provides high quality processing of noisy DWFS data automatically. No ad hoc tuning of the regularization constant is necessary. CLS object estimates are compared with those processed manually. In all cases, the new technique provides images with comparable resolution. In addition, the algorithm is computationally inexpensive, converging to a solution in less than 10 iterations. The remainder of this paper is organized as follows. An image degradation model is presented in Section 2 to include the effects of shot noise and detector read noise. Section 3 reviews the traditional DWFS estimator. In Section 4, we present a CLS algorithm which incorporates a new constraint based on the ensemble average DWFS data. Section 5 presents sample results and Section 6 provides a brief conclusion.
2. Image model The standard model for a linear shift-invariant imaging system can be written as w23x i Ž x , y . s hŽ x , y .) o Ž x , y . ,
Ž1.
where iŽ x, y . is the image, hŽ x, y . is the impulse response or PSF, oŽ x, y . is the object irradiance, and ) denotes two-dimensional convolution. Taking the Fourier transform of both sides of Eq. Ž1. yields I Ž u,Õ . s H Ž u,Õ . O Ž u,Õ . ,
Ž2.
where a capital letter denotes the Fourier transform of an associated lower case quantity and H Ž u,Õ . is the Fourier transform of hŽ x, y ., also known as the OTF. In this paper, the OTF is a random quantity representing the combined effect of atmospheric turbulence and fixed telescope aberrations. Atmospheric turbulence causes a random phase aberration in the telescope pupil. The corresponding OTF is given by a normalized autocorrelation of the generalized pupil function and is a random process due to the presence of the turbulence-induced phase aberrations w23x. The OTF for a diffraction-limited full aperture is a deterministic tapered low-pass filter. Random aberrations in the telescope pupil further degrade performance resulting in a more attenuated OTF, especially at high spatial frequencies w23x. The result is a significant reduction in resolution when compared to the diffraction-limit. For large uncompensated telescopes collecting long exposure images, the resolution is approximately equal to a telescope of diameter r 0 , where r 0 is the Fried parameter which reflects turbulence strength w2x. Since typical values of r 0 range from less than 10 cm to 20–30 cm, resolution is usually limited by atmospheric turbulence not the optical system. To further complicate matters, the image iŽ x, y . is also degraded during the detection process. CCD detectors, which are often used to collect DWFS data, exhibit noise effects due to the random arrival time and location of photoevents. This noise is known as shot or photon noise and is governed by Poisson statistics w2x. CCD detectors also suffer from read noise effects denoted by the quantity n r Ž x, y .. In this paper, we assume the read noise is zero mean, Gaussian distributed, and spatially uncorrelated with uniform variance sr2. Given these assumptions, a detected image dŽ x, y . can be modeled using w3x K
dŽ x , y. s
Ý d Ž x y x k , y y yk . q n r Ž x , y . ,
Ž3.
ks 1
where K is the number of photoevents in the image, Ž x k , y k . is the location of the k th photoevent, and d denotes the Dirac delta function. In this model, K, Ž x k , y k ., and n r Ž x, y . are random variables. Since detector read noise is a form of signal independent noise, n r Ž x, y . is independent of K, Ž x k , y k ., and H Ž u,Õ ..
S.D. Ford et al.r Optics Communications 151 (1998) 93–100
To investigate the modified CLS formulation, a vectormatrix expression is needed for Eq. Ž3.. The new expression includes an additive noise vector n which incorporates both shot and read noise effects based on writing the total noise as nŽ x, y . s dŽ x, y . y iŽ x, y .. Thus, Eq. Ž3. can now be written as d s Ho q n,
Ž4.
where d and o are P length vector versions of the functions dŽ x, y . and oŽ x, y ., respectively. P is the number of pixels in the detector array. The matrix H is a P = P block-circulant matrix representing the shift-invariant PSF. H and o are properly ordered to perform the discrete convolution of Eq. Ž1. w20,21x.
95
4. Modified constrained least-squares formulation Our CLS estimation problem can be stated as follows: ‘‘Given the DWFS-derived estimate ² H˜i) Ž u,Õ . Di Ž u,Õ .:, the ensemble average magnitude-squared OTF estimate ²< H˜i Ž u,Õ .< 2 :, and statistical models for the OTF and noise, find the CLS optimal estimate of the object that caused the data.’’ To accomplish this task, consider finding an object estimate oˆ that minimizes ²< < H˜i Coˆ < < 2 : ,
Ž7.
subject to the constraint < <² H˜iTd i : y ² H˜iTH˜i : oˆ < < 2 s E < < H T n < < 2 4 ,
3. Traditional estimator Raw DWFS data consists of an ensemble of short exposure images and corresponding estimates of the wave-front phase from a WFS. The phase estimates f˜ i Ž x, y . can be used to generate an estimate of the OTF H˜i Ž u,Õ . via a normalized autocorrelation of the generalized pupil function where the subscript i refers to the ith realization. The DWFS estimator suggested by Primot et al. w11x can be written as Oˆ Ž u,Õ . s
² H˜i) Ž u,Õ . Ii Ž u,Õ . : 2
¦ H˜ Ž u,Õ . ;
,
Ž5.
i
where ² : denotes an ensemble average. Eq. Ž5. does not provide acceptable reconstructions due to residual noise in the OTF estimation process. When detector noise is present, Ii Ž u,Õ . is not available and Di Ž u,Õ . must be substituted in Eq. Ž5.. The standard solution is to incorporate a manually selected regularization constant e in the denominator of Eq. Ž5. which gives a new estimator w13,14x Oˆ Ž u,Õ . s
² H˜i) Ž u,Õ . Di Ž u,Õ . : 2
¦ H˜ Ž u,Õ . ;q e
.
Ž6.
i
The regularization constant e is adjusted manually by the user based on the perceived image quality. In the next section, we investigate an algorithm which replaces manual selection of e with an automatic algorithm based on CLS estimation. The algorithm takes advantage of noise reduction through ensemble averaging by directly incorporating ² H˜i) Ž u , Õ . D i Ž u , Õ .: and ²< H˜i Ž u,Õ .< 2 :. Thus, the objective function used to derive this estimator takes on an unfamiliar form when compared to more traditional CLS applications in image processing w19–21x.
Ž8.
where C is a block-circulant constraint matrix which enforces prior knowledge of the true object such as smoothness or support, H˜i is the block-circulant estimated PSF matrix for the ith realization, oˆ is an ordered vector containing the object estimate, and < < < < 2 is the normsquared operator. Notice that this problem statement is different from the form seen in previous CLS image p ro cessin g ap p licatio n s w 1 9 – 2 1 x . H ere, th e ² H˜i) Ž u,Õ . Di Ž u,Õ .: and ²< H˜i Ž u,Õ .< 2 : quantities are incorporated in the constraint function directly. Also, the mean-square error associated with the measurement is now E< < H T n < < 2 4 instead of the norm-squared of the noise E< < n < < 2 4 used previously w19–21x. Even though Eqs. Ž7. and Ž8. are unfamiliar, we will show that the associated Fourier domain filter has the same general form as Eq. Ž6.. The appropriate objective function J Ž oˆ . for the Lagrange minimization can be written as J Ž oˆ . s oˆT C T² H˜iTH˜i : Coˆ q l < <² H˜iTd i : y² H˜iTH˜i : oˆ < < 2 y E < < H T n < < 2 4 ,
4
Ž9.
where l is a Lagrange multiplier w20,21x. We can minimize the objective function in Eq. Ž9. by differentiating J Ž oˆ . with respect to o, ˆ setting the derivative equal to zero, and solving for o. ˆ Applying the appropriate vector-matrix identities w24x to take the derivatives in Eq. Ž9. yields y1
oˆ s ² H˜iT H˜i : ² H˜iT H˜i : q g C T² H˜iT H˜i : C
ž
=² H˜iTH˜i : ² H˜iTd i : ,
/ Ž 10.
where g s 1rl. It is not convenient to evaluate Eq. Ž10. directly since C and H˜i are large matrices for standardsized detector arrays. However, these matrices are block circulant. Block circulant matrices are diagonalized by the discrete Fourier transform allowing the transform domain
S.D. Ford et al.r Optics Communications 151 (1998) 93–100
96
Rewriting the E< < H T n < < 2 4 term in the Fourier domain using Parseval’s theorem w21x yields E < < H T n < < 2 4 sE
½Ý
H ) Ž u,Õ . N Ž u,Õ .
u ,Õ
sE
½Ý ½Ý
2
5
H ) Ž u,Õ . Ž I Ž u,Õ . y D Ž u,Õ . .
2
u ,Õ
sE
5
H ) Ž u,Õ . Ž H Ž u,Õ . O Ž u,Õ . y D Ž u,Õ . .
2
u ,Õ
5
.
Ž 13. Substituting the Fourier transform of the noise model given in Eq. Ž3. into Eq. Ž13. introduces the noise statistics such that E < < H T n < < 2 4 sE Fig. 1. Satellite computer rendering used to test CLS algorithm performance. Negative image shown for clarity.
½
Ý
ž
H ) Ž u,Õ . KH Ž u,Õ . On Ž u,Õ .
u ,Õ
2
K
y
Ý exp wyj2p Ž ux k q Õyk . x y Nr Ž u,Õ . ks 1
equivalent of Eq. Ž10. to provide a simpler scalar form w20,21x. Cancelling ²< H˜i Ž u,Õ .< 2 : terms yields Oˆ Ž u,Õ . s
² H˜i) Ž u,Õ . Di Ž u,Õ . : 2
¦ H˜ Ž u,Õ . ;q g C Ž u,Õ .
2
,
Ž 11.
i
where C Ž u,Õ . is the Fourier transform of the constraint cŽ x, y . associated with the matrix C. For all results presented in Section 5, C is the identity matrix and C Ž u,Õ . s 1 for all spatial frequencies Ž u,Õ .. Thus, Eq. Ž11. is identical to Eq. Ž6.. In this case, our iterative CLS estimation scheme provides automatic selection of the regularization constant. However, Eq. Ž11. will accommodate a general Fourier domain constraint function which is tailored to prior knowledge or a specific application. Instead of tuning g manually, we use a Newton-Raphson technique to find the parameter value which minimizes the objective function J Ž oˆ . given in Eq. Ž9.. The most efficient iteration scheme is on the Lagrange multiplier l not its inverse g w19x. To implement the Newton-Raphson technique, we require the derivative of J Ž oˆ . with respect to l. This derivative is straightforward and given by
E J Ž oˆ . El
s < <² H˜iTd i : y ² H˜iTH˜i : oˆ < < 2 y E < < H T n < < 2 4 .
Ž 12.
The first norm-squared term in Eq. Ž12. is a function of the ensemble average data and the current object estimate. It can be computed conveniently using the Fourier domain quantities ² H˜i) Ž u,Õ . Di Ž u,Õ .: and ²< H˜i Ž u,Õ .< 2 :. However, the second norm-squared term E< < H T n < < 2 4 is not a function of DWFS data and must be derived. Here, we will make use of the noise model given in Section 2.
/5
,
Ž 14.
where the normalized object spectrum OnŽ u,Õ . s O Ž u,Õ .r K has been introduced to the above expression and K is the average number of photoevents per image. Now the linearity property allows us to move the expectation inside the summation over spatial frequencies Ž u,Õ . and evaluate one term of Eq. Ž14.. Standard techniques for evaluating expectations of doubly stochastic Poisson random processes use nested conditional expectations over the random quantities Ž x k , y k ., K, and H Ž u,Õ . w3x. After cancelation of terms, the result for a single spatial frequency is E < H ) Ž u,Õ . N Ž u,Õ . < 2 4 s Ž Kq Psr2 . E < H Ž u,Õ . < 2 4 .
Ž 15. Summing the result in Eq. Ž15. over all spatial frequencies gives the final result E < < H T n < < 2 4 s Ž Kq Psr2 . Ý E < H Ž u,Õ . < 2 4 .
Ž 16 .
u ,Õ
If the OTF statistics are unavailable in Eq. Ž16., the OTF estimate data can be substituted which yields E < < H T n < < 2 4 f Ž Kq Psr2 . Ý ² < H˜i Ž u,Õ . < 2 :.
Ž 17 .
u ,Õ
The result given in Eq. Ž17. is used in Eq. Ž12. to find the current derivative of J Ž oˆ . with respect to the Lagrange multiplier l. This derivative is then used in the NewtonTable 1 Average photoevents per integration time mn
KW
K
4 6 8
3234 518 83
166469 26635 4262
S.D. Ford et al.r Optics Communications 151 (1998) 93–100
97
Fig. 2. CLS algorithm comparison with manual parameter selection, r 0 s 20 cm, mn s q4, sr s 15 electrons per pixel. Ža. Single short exposure image. Žb. CLS algorithm estimate, g s 0.0015, 4 iterations. Žc., Žd., Že., Žf. Manual parameter selection, e s 0.1, 0.01, 0.0001, and 0.00001, respectively.
98
S.D. Ford et al.r Optics Communications 151 (1998) 93–100
Fig. 3. CLS algorithm output versus object brightness, r 0 s 10 cm, sr s 15 electrons per pixel. Ža. Detected data, mn s q4. Žb. CLS estimate, mn s q4, g s 0.00027, 4 iterations. Žc. Detected data, mn s q6. Žd. CLS estimate, mn s q6, g s 0.0093, 5 iterations. Že. Detected data, mn s q8. Žd. CLS estimate, mn s q8, g s 0.8180, 7 iterations.
S.D. Ford et al.r Optics Communications 151 (1998) 93–100
Raphson iteration to update the value of l and generate an object estimate using Eq. Ž10. or its Fourier domain equivalent Eq. Ž11.. The iteration continues until the object estimate meets a pre-determined stopping criterion. A variety of criteria can be used to stop the Newton-Raphson iteration w22x. Since our goal is to minimize the objective function given in Eq. Ž9., we choose to stop the iteration when J Ž oˆ . is sufficiently small. All CLS images were processed until J Ž oˆ . F 0.001. In the next section, we show the results of using this technique on simulated DWFS data.
5. Sample results The sample results shown in this section are based on a computer rendering of a typical satellite as shown in Fig. 1. The detector array was 256 = 256 elements for a total of P s 65536 pixels. The satellite was 10 m in length and in low earth orbit at a range of 500 km. The diameter of the telescope was 1 m with both imaging and wavefront sensor wavelengths set at ls 500 nm. The telescope fractional bandwidth was "5%. The object was assumed to have the same spectral signature as the sun. The WFS subaperture size was 10 cm for a total of 60 subapertures within the entrance pupil. Atmospheric turbulence and detector noise effects were modeled using an existing computer simulation w3x. For each data realization, the simulation created a random phase screen with the appropriate turbulence statistics w25x, calculated the true OTF, and formed a detected image realization d i Ž x, y .. At the same time, a WFS model was used to generate a phase estimate f˜ i Ž x, y . which was then used to compute the estimated OTF. Finally, the required quantities H˜i) Ž u,Õ . Di Ž u,Õ . and < H˜i Ž u,Õ .< 2 were accumulated and the process repeated 150 times to generate the ensemble averages. The simulation incorporates an intensity splitter which sends 40% of the photons to the image plane and 60% to the WFS. Integration times of 10 ms were used for both WFS and imaging system. We modeled the object brightness associated with visual magnitudes mn s q4, q6, and q8. The resulting photoevents per integration time for these cases are presented in Table 1, where KW is the average number of photoevents across an individual WFS subaperture. Fig. 2 displays the detected short exposure image and CLS algorithm outputs associated with excellent seeing conditions Ž r 0 s 20 cm., a relatively bright satellite object Ž mn s q4., and detector read noise representing a high quality CCD detector Ž sr s 15 electrons per pixel.. A single short exposure image is given in Ža. to illustrate shot noise and detector read noise effects. Clearly, this image is degraded by the detection process. After optimal CLS processing of the DWFS data as shown in Žb., image resolution is greatly enhanced. In this case, the CLS
99
algorithm automatically selected the regularization parameter g s 0.0015 in 4 iterations. In Žc., Žd., Že., and Žf. the regularization parameter e associated with the traditional DWFS estimator given in Eq. Ž6. was selected manually. For images Žc. and Žd., e ) g . For images Že. and Žf., e - g . Notice that the quality of the CLS image in Žb. is as good or better than the manual images. Thus, Fig. 2 illustrates that the CLS algorithm selected a reasonable value for g in this case. Not only does shot noise degrade the short exposure data, it also restricts WFS accuracy. Without a sufficiently accurate OTF estimate, DWFS performance is severely degraded. To illustrate this limitation while applying the CLS algorithm, consider Fig. 3. Here, r 0 s 10 cm and detector read noise remains unchanged from the previous cases. Images Ža., Žc., and Že. show single short exposure realizations for object brightness cases mn s q4,q 6, and q8, respectively. In Žc. and Že., light level is so low, detector read noise is clearly dominant upon visual inspection of the images. Images Žb., Žd., and Žf. give the corresponding CLS estimates where g s 0.00027, 0.0093, and 0.8180, respectively. As object brightness decreases, the output images are more blurred. This effect is consistent with deconvolution using a poor quality estimate of the OTF and noise effects driving the CLS regularization parameter selection to a value near unity.
6. Conclusion We investigated a CLS estimator that incorporates noisy DWFS data, noise statistics, and OTF statistics. This estimator selects the regularization parameter automatically. No ad hoc tuning is necessary to reduce high spatial frequency noise effects in the DWFS image. The CLS estimator uses a Newton-Raphson iteration to select a Lagrange multiplier which minimizes an objective function. The objective function uses ensemble average data directly which aids in noise suppression. Our simulation results show that the new algorithm produces DWFS images comparable in quality to manual regularization with minimal computational expense. Finally, the CLS estimator was derived to incorporate a general Fourier domain constraint. Thus, other constraint functions can be used based on the specifics of the application.
Acknowledgements This work was supported by the Air Force Office of Scientific Research, Bolling AFB, Maryland and the Air Force Maui Optical StationrMaui High Performance Computer Center Research and Development Consortium in Kihei, Maui, Hawaii.
100
S.D. Ford et al.r Optics Communications 151 (1998) 93–100
References w1x F. Roddier, The Effects of Atmospheric Turbulence in Optical Astronomy, in: E. Wolf ŽEd.., Progress in Optics-XIX, North-Holland, New York, 1981. w2x J.W. Goodman, Statistical Optics, Wiley, New York, 1985. w3x M.C. Roggemann, B.M. Welsh, Imaging Through Turbulence, CRC Press, Boca Raton, Florida, 1996. w4x J.H. Hardy, Proc. IEEE 66 Ž1978. 651. w5x A. Labeyrie, Astron. Astrophys. 6 Ž1970. 85. w6x K.T. Knox, B.J. Thompson, Astrophys. J. 193 Ž1974. 45. w7x A.W. Lohmann, G. Weigelt, B. Wirnitzer, Appl. Optics 22 Ž1983. 4028. w8x M.C. Roggemann, D.W. Tyler, M.F. Bilmont, Appl. Optics 31 Ž1992. 7429. w9x M.C. Roggemann, C.A. Stoudt, B.M. Welsh, Opt. Eng. 33 Ž1994. 3254. w10x J. Fontanella, J. Opt. ŽFrance. 16 Ž1985. 257. w11x J. Primot, G. Rousset, J.C. Fontanella, J. Opt. Soc. Am. A 7 Ž1990. 1598. w12x D.L. Fried, in: Proc. 828 SPIE Conference on Digital Image Recovery and Synthesis, 1987, pp. 127–133. w13x J.D. Gonglewski, D.G. Voelz, J.S. Fender, D.C. Dayton, B.K. Spielbusch, R.E. Pierson, Appl. Optics 29 Ž1990. 4527.
w14x B.M. Welsh, R.N. VonNiederhausern, Appl. Optics 32 Ž1993. 5071. w15x M.C. Roggemann, B.M. Welsh, J. Devey, Appl. Optics 33 Ž1994. 5754. w16x T.J. Schulz, in: Proc. 2029 SPIE Conference on Digital Image Recovery and Synthesis II, 1993, pp. 3111–320. w17x M.C. Roggemann, B.M. Welsh, Appl. Optics 33 Ž1994. 5400. w18x B.M. Welsh, M.C. Roggemann, Appl. Optics 34 Ž1995. 2111. w19x M.C. Roggemann, D.W. Tyler, Appl. Optics 36 Ž1997. 2360. w20x B.R. Hunt, IEEE Trans. Comput. C 22 Ž1973. 805. w21x R.C. Gonzales, R.E. Woods, Digital Image Processing, Addison-Wesley, Reading, Massachusetts, 1992. w22x E. Isaacson, H.B. Keller, Analysis of Numerical Methods, Dover Publications, New York, 1994. w23x J.W. Goodman, Introduction to Fourier Optics, McGraw-Hill, New York, 1968. w24x J.L. Melsa, D.L. Cohn, Decision and Estimation Theory, McGraw-Hill, New York, 1978. w25x B.M. Welsh, A Fourier Series based atmospheric phase screen generator for simulating nonisoplanatic geometries and temporal evolution, submitted to Appl. Optics, February 1997, unpublished.