Image reconstruction with analytical point spread functions
A. Asensio Ramos^{1,2}  A. López Ariste^{3}
1  Instituto de Astrofísica de Canarias, 38205 La Laguna, Tenerife, Spain
2 
Departamento de Astrofísica, Universidad de La Laguna, 38205 La Laguna, Tenerife, Spain
3 
THEMIS, CNRS UPS 853, c/vía Láctea s/n, 38200 La Laguna,
Tenerife, Spain
Received 24 February 2010 / Accepted 19 April 2010
Abstract
Context. The image degradation produced by atmospheric
turbulence and optical aberrations is usually alleviated using
postfacto image reconstruction techniques, even when observing with
adaptive optics systems.
Aims. These techniques rely on the development of the wavefront
using Zernike functions and the nonlinear optimization of a certain
metric. The resulting optimization procedure is computationally heavy.
Our aim is to alleviate this computational burden.
Methods. We generalize the extended ZernikeNijboer theory to
carry out the analytical integration of the Fresnel integral and
present a natural basis set for the development of the point spread
function when the wavefront is described using Zernike functions.
Results. We present a linear expansion of the point spread
function in terms of analytic functions, which, in addition, takes
defocusing into account in a natural way. This expansion is used to
develop a very fast phasediversity reconstruction technique, which is
demonstrated in terms of some applications.
Conclusions. We propose that the linear expansion of the point
spread function can be applied to accelerate other reconstruction
techniques in use that are based on blind deconvolution.
Key words: techniques: image processing  methods: analytical  methods: numerical  telescopes  atmospheric effects
1 Introduction
Atmospheric turbulence degrades astronomical images by introducing aberrations in the wavefront. In the past few years, functional adaptive optics systems have been developed with the aim of partially cancelling these aberrations and generating a diffractionlimited image of any astronomical object (e.g., Beckers 1993). The development of these systems has been particularly challenging for solar telescopes because the reference object (the solar surface) has spatial structure. In spite of the great success of adaptive optics systems (Scharmer et al. 2000; Rimmele 2000), the limited number of degrees of freedom of the deformable mirror and the speed at which it has to work, especially in solar observations, impedes a full correction of the wavefront. Therefore, it has become customary in solar imaging to apply postfacto reconstruction techniques to the observed images to improve the wavefront correction. The advantage of these techniques is that the time limitation is much less restrictive because they do not need to work in realtime. Hence, one can use lots of computation power to reconstruct the images. One of the methods first considered for the reconstruction of solar images was based on speckle techniques (von der Lühe 1994,1993), which produced images of extraordinary quality. However, it is affected by two fundamental problems. On the one hand, the number of images that one needs to acquire is large (although comparable to the number of images needed in complex blind deconvolution algorithms) because the missing information about the wave phases is statistically estimated from a succession of many short exposures. On the other hand, it relies on a theoretical description of the effect of the atmospheric turbulence and adaptive optics on the image (Wöger & von der Lühe 2007), which might be inaccurate. However, studies have shown that almost realtime processing (Denker et al. 2001) and photometric precision (Wöger et al. 2008) can be achieved routinely.A complementary approach to speckle reconstruction appeared with the systematic application of techniques based on phase diversity (e.g., Löfdahl & Scharmer 1994a; Löfdahl et al. 1998), developed some years before by Gonsalves & Chidlaw (1979) and later extended by Paxman et al. (1992). These techniques rely on a maximum likelihood simultaneous estimation of the wavefront and the original image (before being affected by atmospheric+telescope aberrations) using a reference image and a second one that contains exactly the same aberrations as the reference one plus a known artificiallyinduced aberration. Because of the mechanical simplicity, a defocus typically constitutes the added aberration. The main drawback of the phase diversity technique is the heavy computational effort needed to compute the maximum likelihood solution. Many Fourier transforms have to be used to estimate the point spread function (PSF) from the wavefront at the pupil plane.
Even more computationally expensive is the blind deconvolution technique developed by van Noort et al. (2005). This technique, termed multiobject multiframe blind deconvolution (MOMFBD), is able to reconstruct estimates of the object and the wavefront from the observation of different objects (same location on the Sun observed at different wavelengths) with some constraints, which are essentially used to augment the amount of information introduced into the likelihood function. As a consequence, the estimation of the wavefront and the original image is less dominated by noise and hence more stable. In close analogy with phase diversity, this enhancement in stability comes at the price of a higher computational burden since many Fourier transforms have to be computed for the maximization of the likelihood. In conclusion, all these techniques depend on timeconsuming computations. As modern solar telescopes increase their flux of high quality data, observers are hampered by the bottleneck of reconstruction algorithms.
Nijboer (1942) made it possible to estimate the PSF of an aberrated optical device in terms of the generalized complex pupil function. The resulting integral expressions were very complicated and for a long time remained impossible to evaluate. What is presently known as the extended NijboerZernike theory (ENZ, Janssen 2002; Braat et al. 2002) presents analytical expressions for the complex field at the focal volume that depend on the aberration described by the coefficients of the Zernike expansion of the wavefront. The expressions are closed in the sense that they depend only on computable functions and the coefficients of the Zernike expansion.
Our purpose in this paper is to extend the ENZ theory to arbitrary aberrations exploiting the mathematical properties of the Zernike functions. As a consequence, we find a general expression for the PSF in the focal volume as a linear combination of given basis functions. Starting from this expression, we propose to use it in a phase diversity reconstruction algorithm. Because of the analytical character of the PSF, we are able to avoid the calculation of Fourier transforms during the iterative process, providing an important gain in computational time. This is our fundamental result: that by applying this formalism we can drastically diminish the computational burden of reconstruction algorithms and therefore the time needed to perform the reconstruction. We illustrate this using phase diversity. In principle, any reconstruction algorithm can, however, benefit from our formalism, since we only have changed the way in which the PSF is described.
The outline of the paper is the following. Section 2 describes the wavefront expansion in terms of Zernike functions and how the field distribution can be derived analytically for general aberrations. The expansion of the point spread function is presented in Sect. 3 with an statistical analysis of that expansion for atmospheric Kolmogorov turbulence in Sect. 4. Section 5 presents the application of the formalism to the rapid restoration of images using phasediversity. Finally, our conclusions are presented in Sect. 6.
2 Wavefront expansion
We define
to represent the wavefront of a point source beam in a pupil
plane, where
and
are the polar coordinates of the unit
disk in that plane. The field
distribution at an image plane can be written in the Fresnel framework as (Born & Wolf 1980)
which is, in general, a complex quantity. In this expression, we define a coordinate system over the image space, where are angular coordinates in the focal plane and f is the distance to that focal plane. We define the focal plane as the plane where the optical system would form the image of the object in the absence of any aberration. We implicitly assume an optical system with numerical aperture , without loss of generality for solar instruments. The coordinates can be related to the geometrical coordinates (X,Y,Z) of the optical system by means of the following relationships
x  =  
r  =  (2) 
where is the wavelength of the light.
2.1 Wavefront expansion
To solve the integral in Eq. (1), we have to provide a useful
description of the wavefront
.
In general, the wavefront is
also a complex quantity that is usually represented in polar form: the modulation in amplitude
is represented by
,
while the phase modulation is determined by
.
If the pupil plane is fully transmissive (or perfectly reflective), we can assume
so that
(3) 
This is the case for a clear aperture telescope where the wavefront is modified by atmospheric turbulence. For the sake of simplicity, we assume this to be the case, since it produces simpler mathematical expressions that can be related more transparently to previous results. An annular pupil would only require a change in the integration limits of Eq. (1), thus modifying accordingly the final result.
It is customary (and for many reasons interesting) to express the phase modulation
of the wavefront in terms of the Zernike functions (Noll 1976).
The Zernike functions are labeled by two quantum numbers: a principal number nand an azimuthal number m fulfilling the conditions that
and that nm is even. Their definition
is
=  
=  
=  (4) 
where is a Zernike polynomial (e.g., Born & Wolf 1980). To simplify mathematical manipulations, Noll (1976) introduced a useful single ordering index j, which has perfect correspondence with pairs so that higher values of j represent higher aberrations with smaller probability amplitudes in a Kolmogorov turbulent atmosphere. We use the two notations interchangeably. The reader should be aware that any pair (n, m) appearing explicitly in an expression has to be in agreement with the corresponding Noll index j. Since Zernike functions constitute a family of orthogonal functions in the unit circle, any phase aberration can be written as a linear expansion
In a turbulent atmosphere of Kolmogorov type, the aberrations are characterized by coefficients a_{j} that follow a multivariate Gaussian statistical distributions whose covariance matrix can be obtained analytically (Noll 1976). It is interesting to point out that the covariance matrix of the coefficients is not diagonal, although it is possible to find modifications of the Zernike functions that diagonalize it. Returning to the expression of the exponential of the phase aberration, we expand it as a power series:
(6) 
The radius of convergence of this expansion is infinite and that does not translate into an easy cutoff of the series. But using the decomposition in terms of Zernike functions of the phase aberration, we end up with
For small aberrations (), the previous series can be safely cut at first order, something that was successfully exploited by Janssen (2002) and Braat et al. (2002) to empirically estimate point spread functions of microscopetype optical systems. This approximation may be useful if one is interested in describing aberrations introduced by well adjusted optical systems, but it is hardly acceptable for aberrations introduced by atmospheric turbulence. Because of the complexity of the resulting expression, there has been no interest in the past in this approach as a means of estimating the PSF of atmospheric seeing. The reason is that the series expansion of Eq. (7) generates products Z_{j} Z_{k} at second order and highorder products at higher orders. In the absence of any clear rule for manipulating these products, the power expansion was rendered useless already at second order.
The mathematical result of Mathar (2009)^{} has radically simplified the
approach to that expansion and we use those results here for the first time to pursue the
analytical integration of the PSF of any aberration introduced by
atmospheric turbulence. Mathar (2009) constructively demonstrated that the product of two Zernike functions
can be written as a linear combination of Zernike functions (see Appendix A), so that
where the coefficients d_{k} can be obtained by direct application of the relations shown in Appendix A. The product expansion of Mathar (2009) can be applied recursively without difficulties, resulting in
(9) 
where the c_{k} coefficients can be calculated from the d_{k} coefficients of Eq. (8). As a consequence, the expansion of Eq. (7) can be written confidently as
where the coefficients in general are complex. In other words, if the phase aberration can be expanded in Zernike functions, the wavefront can also be expanded efficiently in these functions with complex coefficients that can be inferred from the coefficients a_{j} by a repetitive application of the coupling relation of Eq. (A.2) and the coupling coefficient of Eq. (A.3). A direct consequence of Eqs. (10) and (7) is that, for weak turbulence, the imaginary parts of the coefficients are approximately equal to the Zernike coefficients of the expansion of the wavefront
(11) 
and strictly equal at first order. The weak turbulence condition for that approximation often applies to our astronomical applications as demonstrated statistically in Sect. 4. By means of this approximation, any method that allows us to estimate the value of the coefficients from the PSF (such as the phase diversity method shown below) provides us directly with an estimation of the Zernike coefficients of the wavefront, at least to first order.
2.2 The field distribution
Having written the pupil function
as a linear combination
of Zernike functions with complex coefficients ,
the
field distribution of Eq. (1) simplifies to a sum of integrals of
the form
In the usual mathematical development of the ENZ theory, we decouple the two integrals above into radial and azimuthal parts. The azimuthal part can be simplified using the JacobiAnger identity (Abramowitz & Stegun 1972)
(13) 
where the are Bessel functions of the first kind. Making use twice of the previous integral identity in Eq. (12) and reformulating the azimuthal part of the Zernike functions into a single (that contains both the and functions into its real and complex parts for an unsigned m index), we end up with
(14) 
The function V_{n}^{m}(r,f) represents the following (complex) radial integral
which fulfills the condition
(16) 
If we are interested in the field distribution at the plane where the image would form in the absence of any optical aberration (f=0), the radial integral becomes real and has a straightforward solution:
(17) 
Prefiguring its application in phase diversity, the integral of Eq. (15) can be numerically calculated for arbitrary values of f using the rapidly convergent series expansion developed by Braat et al. (2002).
The final expression for the field distribution is given by
We note that, when no aberrations are present, the infocus field distribution (at f=0) is proportional to
(19) 
which is the function giving rise to the well known Airy spot for a circular nonaberrated aperture. At , the field distribution is proportional to V_{0}^{0}(r,f), which, using the rapidly convergent series of Braat et al. (2002), can be written exactly as
(20) 
To summarize, the product expansion of Mathar (2009) has allowed us to write the field distribution in the image space as a linear combination of functions V_{n}^{m}(r,f) defined above, and of which the (infocus or outoffocus) Airy spot is just the first contribution. This can be considered to be a generalization of the results presented by Janssen (2002) and Braat et al. (2002) for the case of small aberrations.
3 The point spread function
The point spread function or intensity distribution is obtained from the complex field distribution as
(21) 
For the sake of simplicity in the presentation, we separate the analysis of the infocus and the outoffocus PSFs.
3.1 Infocus PSF
Using Eq. (18) with the first term included in the summation being j=0, the PSF results in the double summation(22) 
where the (n,m) pair is associated with the Noll index j, while the (n',m') pair is attached to the Noll index k. The previous double summation can be separated into two summations: one containing terms for which j=k and one containing terms . After tedious but straightforward algebra, the final expression for the point spread function is
The functions appearing in the linear expansion can be calculated as
We note that the azimuthal and radial contributions are separate. The angular dependence is modulated by the functions
The coefficients C_{jk} and D_{jk} in the linear expansion of the PSF are related to the real and imaginary parts of products of coefficients
C_{jk}  =  
D_{jk}  =  (26) 
where . As an example, Fig. 1 shows the unique basis functions up to Noll index j=7. The remaining functions with different combinations of indices can be related to these by trivial relations such as those shown in Appendix B. We note that the terms j=k of the summation contain only the contribution of the C_{jk} coefficients and are azimuthindependent, since
=  
=  0.  (27) 
Because nm is even in the definition of the Zernike functions, it can be verified that only the terms with j=k in Eq. (23) contribute to the total area of the PSF, so that
where the value of n is associated with the Noll index of summation. As a result, the Strehl ratio defined as the ratio of the peak intensity at the origin of the PSF to the equivalent Airy function simplifies to
(29) 
Figure 1: Example of unique basis functions for representing the PSF, including terms up to the Noll index j=7. These functions are calculated for a telescope of 100 cm diameter at a wavelength of 5250 Å, a typical situation in presentday telescopes. The spatial dimensions are in units of pixels, that we choose to be of size 0.055 , also typical of present observing conditions. The upper panel shows the results for the focused image while the lower panel shows the results for the defocused image, where the defocus is identical to that shown in Sect. 5. 

Open with DEXTER 
3.2 Outoffocus PSF
Following the same approach as in the previous section, we can also write an expression for the PSF at any arbitrary focal position around the focal point. Since in this case the V_{n}^{m}(r,f) functions in general are complex, the basis functions are slightly more elaboratewhere
The auxiliary functions appearing in the previous formulae are given by
(32) 
When the previous expressions are evaluated at the focus (f=0), because the V_{n}^{m} functions are real, we find of course that and , thus recovering the expressions of Eq. (24). Likewise, for the terms of the summation in Eq. (30) for which j=k, we have
=  
=  0.  (33) 
Figure 2: Histograms of the distribution of the wavefront expansion coefficients for an atmosphere with Kolmogorovtype turbulence (upper row). The real and imaginary part of the corresponding expansion coefficients of the field distribution are shown in the middle and lower rows, respectively. The simulations are carried out for a telescope of D=90 cm and a Fried parameter r_{0}=7 cm. 

Open with DEXTER 
4 Statistical properties of the coefficients
In a Kolmogorov turbulent atmosphere, the coefficients of the Zernike expansion of the wavefront follow a multivariate Gaussian distribution with a nondiagonal covariance matrix determined by Noll (1976), Wang & Markey (1978), and Roddier (1990). The amplitude of the elements of the covariance matrix is determined by the Fried parameter r_{0}. As a function of this parameter, we simulate random samples from this multivariate Gaussian and use the expressions of Sect. 2 and the formulae developed by Mathar (2009) presented in Appendix A to calculate the equivalent coefficients. The corresponding statistical analysis is presented in Fig. 2 for aberrations of order higher than two, excluding tiptilt. The upper row shows that the coefficients are Gaussian distributed, with a variance depending on the order of the aberration. The middle and lower rows show the histograms for the real and imaginary part of the coefficients, respectively.
The distribution of the imaginary part is close to Gaussian and very similar to that of the wavefront expansion coefficients. As already mentioned, this is a consequence of, to first order, the imaginary part of the coefficients being equal to the coefficients. Higherorder corrections to the imaginary parts occur only at odd orders, so that the first correction already takes place at thirdorder, which often induces a relatively small change for turbulence that is not too strong. The Pearson linear correlation coefficient between and is larger than 0.95 for practically all orders, indicating that represents a good estimation of the wavefront expansion coefficients. The real parts present distributions with a smaller variance and, in some cases, a large kurtosis. This is because the real part of the coefficients is modified only at secondorder often representing a small correction.
5 Phase diversity
A fundamental advantage of the analytical expression for the PSF just presented is that it allows us to introduce a defocus in the PSF, while keeping unchanged the coefficients that characterize the aberrations. As a consequence, it is natural to propose such a functional form when performing image reconstruction first by means of phase diversity techniques (Paxman et al. 1992; Gonsalves & Chidlaw 1979; Löfdahl & Scharmer 1994b). To decide on naming conventions, we first summarize the basics of the standard approach to phase diversity.
5.1 Standard approach
The idea behind phase diversity is to simultaneously acquire a focused image (D_{0}) and another one with exactly the same aberrations plus a known defocus (D_{k}). In the standard image formation paradigm, the observed image after degradation by the atmosphere and the optical devices of the telescope can be written asD_{0}  =  F * P_{0} + N_{0}  
D_{k}  =  F * P_{k} + N_{k},  (34) 
where P_{0} and P_{k} are the PSFs of the atmosphere+telescope in the focused and defocused images, respectively. The quantities N_{0} and N_{k} are noise contributions in the image formation, each one characterized with a variance and , respectively. Since both images are obtained simultaneously, the underlying object F is the same in both cases.
In the standard phase diversity method (see, e.g., Paxman et al. 1992; Gonsalves & Chidlaw 1979; Löfdahl & Scharmer 1994b), the phase aberration
in the focused image is expanded in a Zernike basis according to Eq. (5)^{}. The PSFs are obtained from this phase aberration by calculating
P_{0}  =  
P_{k}  =  (35) 
where the defocus is introduced by adding a term of the form and stands for the inverse Fourier transform. The selection of the defocusing distance is essential for obtaining an efficient and successful phase diversity method. It is usually chosen to represent one wave retardance at the wavelength of the observation.
The coefficients of the expansion in Zernike functions of the phase aberration are obtained by
calculating the coefficients
of the expansion of Eq. (5) that
minimize the error metric (Paxman et al. 1992)
where , , , and are the Fourier transforms of D_{0}, D_{k}, P_{0} and P_{k}, respectively. The summation is carried out over the full Fourier domain defined by the frequencies u and v, while . We note that the dependence of the metric on the coefficients is highly nonlinear because of the exponential function, the Fourier transforms, and the modulus operation that appears in the definition of the PSF. We note also that each evaluation of the metric and its gradient with respect to the a_{j} coefficients involves several Fourier transforms, which are of order . As a final step, once the PSF (and the corresponding Fourier transform) are known, the reconstructed image is estimated as the inverse Fourier transform of (Paxman et al. 1992)
The maximumlikelihood solution of the metric of Gonsalves & Chidlaw (1979) is not affected by the presence of additive Gaussian noise but it strongly affects the deconvolution process given by Eq. (37) since noise is unlimitedly amplified. We follow Löfdahl & Scharmer (1994b) and apply a filter in the Fourier domain to the Fourier transform of the observed images with the aim of filtering a large part of the noise. The same procedures applied by Löfdahl & Scharmer (1994b) to filter the remaining noise peaks at high frequency are also applied.
5.2 The new approach
The crucial advantage of the phase diversity algorithm that we propose is that the PSF of the focused and the defocused images can be written as linear combination of known functions. Because of that, the previous error metric can be minimized without performing a single Fourier transform during the iterative scheme. After Fouriertransforming the observed images and the basis functions X and Y, it is a matter of finding the values of the complex quantities that minimize the metric L. Another advantage is that the nonlinear dependence of the metric function on the coefficients is less pronounced, thus helping the optimization routine to find the minimum.
To find the values of the coefficients
that minimize the metric L, we
use a scaled conjugate gradient algorithm (Andrei 2008) that performs extremely well
and rapidly for complex problems. To this end, one needs the gradient of L with respect to the coefficients. After some tedious but straightforward algebra, the gradient can be found to be
=  
(38) 
where and
=  
=  
=  
=  (39) 
The superindices r and i indicate the real and imaginary part of the coefficient , respectively. One should be aware that the values of n and m in the previous equation are those associated with the Noll index l, while n' and m' are those associated with the index k of the summation. Finally, the functions and are the Fourier transforms of the functions defined in Eq. (31).
Figure 3: Reconstruction example at 395 nm. The lower image shows the infocus image while the right image presents the reconstructed image using only one pair of infocus and defocused images. The contrast of the original image is 9.6% while the reconstructed one presents a contrast of 11.9%. Data was acquired with the SST telescope. 

Open with DEXTER 
The computer code is written in Fortran 90 and runs in parallel environments using the Message Passing Interface library. Prior to the phase diversity restoration procedure, there is an initial phase in which the basis functions of Eqs. (24) and (31) are precalculated and their Fourier transforms are saved for later use. Since this has to be run only once for each telescope and size of the isoplanatic patch, it can be done with the required precision at virtually no extra computing time. In a second phase, in which the actual phase diversity restoration process is carried out, these Fourier transforms of the basis functions are quickly read and used in all subsequent calculations. This greatly accelerates the computing time and only two Fourier transforms are needed per isoplanatic patch. The parallel implementation scales very well with the number of processors. The fieldofview is divided into isoplanatic patches of a given size and the reconstruction of each patch is carried out by one of the available processors. In principle, if processors were at hand, the computing time would be reduced by a factor because the contribution of all input/output tasks would be negligible with respect to the true iterative process for minimizing the phase diversity merit function.
5.3 Illustrative examples
We present two representative examples of the success of the implementation of our analytical representation of the PSF for carrying out phasediversity image reconstruction. We emphasize that our purpose is not to present amazing reconstructions but to demonstrate that reconstruction can be performed to the same level of quality as the standard procedures only much faster.
5.3.1 Blue images
In our first example, the quiet Sun internetwork was observed with the Swedish 1m Solar Telescope (SST; Scharmer et al. 2002) on 2007 September 29 during a period of very good seeing (see Bonet et al. 2008, for a description of the observations). The pixel size is 0.034'' and the observations are carried out in the continuum close to the Ca H line at 3953 Å, providing a theoretical diffraction limit of 0.1''. All these parameters were chosen to enhance granulation contrast. With an isoplanatic patch of pixels ( ), an image of pixels contains approximately 60 isoplanatic patches. Our phase diversity restoration algorithm takes around 5 s per patch when including 20 terms in the summation of Eq. (10) that describes the wavefront. Therefore, the total computing time is 5 min. When operating on a standard presentday desktop computer with 4 processors, every image can be fully restored on a timescale of 1 min. Figure 3 presents the results of the reconstruction (right panel) and the original focused image (left panel). The zoomed region illustrates the spatial frequencies enhanced by the phasediversity algorithm.5.3.2 Red images
In this second example, we analyze observations of an active region carried out during 2002 May 15 with the THEMIS telescope (Rayrole & Mein 1993) at a wavelength of 850 nm. The pixel size is 0.075'' and, due to the long wavelength, the diffraction limit is 0.24'', limiting the contrast of the granulation, both because of the wavelength and the intrinsically reduced contrast in the infrared. This dataset was also used by Criscuoli et al. (2005) to present the first results of phasediversity applied to THEMIS data. The images were restored with patches of pixels, equivalent to . Patches were finally mosaicked to build the final image shown in Fig. 4, where only the part inside the rectangle has been reconstructed, using 20 terms in the summation of Eq. (10) that describes the wavefront. Visual inspection clearly indicates an improvement in the image quality. If desired, higher quality visual and quantitative results can be obtained by adding multiimage information to the reconstruction algorithm (Criscuoli et al. 2005).Figure 4: Example of the phasediversity reconstruction at 850 nm. We represent the focused image (left) and the reconstructed subimage indicated with a box (right). Images were acquired with THEMIS. 

Open with DEXTER 
Computing time is what makes the difference, since each patch, as in the previous case, took only 5 s on a standard desktop computer. With these requirements, image reconstruction can be easily implemented as an online task in the data pipelines of the telescopes themselves.
6 Conclusions
The use of the extended NijboerZernike theory and some mathematical results about the multiplication of Zernike polynomials has allowed us to rewrite the image formation integrals with an aberrated wavefront in terms of linear combinations of analytic functions. Generalizing the usual ENZ theory, we are able to do so a priori independently of the amplitude of the aberrations. With this mathematical tool, we have been able to rewrite the techniques of postfacto image reconstruction taking advantage primarily of the Fourier transforms of those analytic functions being able to be precomputed once and for all. Building different PSFs in the minimization process at the core of these reconstruction algorithms does not require the recomputation of any other Fourier transform but just the modification of the scalar, although complex, coefficients multiplying the functions.
The gain in speed is enormous, as we have demonstrated for two observations prepared for phase diversity. Image reconstruction with a quality identical to the standard algorithms is performed in a matter of seconds per patch, and one to five minutes for widefield images, using standard desktop computers. If image reconstruction were a must for both present and future solar observations, the bottleneck produced by computationally expensive and timeconsuming algorithms would almost certainly be a showstopper for future instruments. The present improvement through the use of analytical PSFs solves the core of that problem and suggests that image reconstruction can even be implemented as a routine online procedure in the data pipelines of the telescope instruments themselves.
Since we have only modified the description of the PSF with respect to previous approaches, all common reconstruction schemes (with any desired complexity), and not just phase diversity, can be extended to use our formalism. Among them, we have identified methods such as multiimage phase diversity (e.g., Paxman et al. 1992) or multiobject multiframe blind deconvolution (van Noort et al. 2005) that use many images to obtain information about different frequencies that have been destroyed by the presence of noise. In essence, all these schemes require writing an error metric (equivalently, a likelihood function) such as Eq. (36) that takes into account the presence of the additional information. The larger amount of information helps us to regularize the problem and reduces the influence of noise. All of the schemes are penalized by the computing time involved and all can potentially profit from the analytical approach we present to reduce the computational burden. The final goal is not to decelerate the flow of images from instruments simply because of out lack of capability to handle the image reconstruction problem.
Beyond those known techniques, we also point out that it is possible to introduce regularization by using a fully Bayesian approach in which prior information is introduced into the problem, eliminating a priori the problem caused by noise in the minimization of the metric.
AcknowledgementsImages from the Swedish Solar Tower (SST) were provided by J. A. Bonet (IAC). Images from THEMIS were provided by D. del Moro (Univ. Roma Tor Vergata). We are very thankful to both of them for providing the data necessary to illustrate this work. Financial support by the Spanish Ministry of Science and Innovation through project AYA200763881 is gratefully acknowledged.
Appendix A: Product of Zernike functions
The breakthrough brought over by Mathar (2009) has been to point out that, in the
product of two Zernike functions, only the following three situations could
arise:
=  
=  
=  (A.1) 
In other words, the result is always the sum of trigonometric functions depending on two azimuthal numbers . This is indeed the seed of two new Zernike functions with if we are able to write the product of the radial polynomials with identical m_{3} azimuthal value. That is, we seek a product expansion for the radial polynomials of the form
Every term of the series above has the same azimuthal number m_{3}, coincident with the m_{3} of the trigonometric product and can therefore be combined with it to build a Zernike function Z_{n3}^{m3}. Mathar (2009) demonstrates that this is indeed feasible and provides the value of the g_{n1,m1,n2,m2,n3,m3} coefficient in the product expansion as
where a=(nm)/2. These coefficients play the same role as the ClebshGordan recoupling coefficients in determining the product of two spherical harmonics.
Appendix B: Some properties of and and analytical Strehl ratio
After interchanging the (n,m) and (n',m') pair of indices, the following expressions hold:
=  
=  
=  
=  (B.1) 
As a consequence, the functions X_{n,n'}^{m,m'} and Y_{n,n'}^{m,m'} fulfill the following properties:
=  
=  (B.2) 
For the Strehl ratio, it is possible to demonstrate that only the Airy part of the PSF (either infocus or outoffocus) contributes to the central peak of the PSF. The reason is that, evaluating Eq. (15) at r=0gives
(B.3) 
so that
=  
=  0,  (B.4) 
which gives, for the f=0, a contribution of 1/4 at the center of the PSF.
References
 Abramowitz, M., & Stegun, I. A. 1972, Handbook of Mathematical Functions (New York: Dover) [Google Scholar]
 Andrei, N. 2008, Optimization, 57, 549 [CrossRef] [Google Scholar]
 Beckers, J. M. 1993, ARA&A, 31, 13 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Bonet, J. A., Márquez, I., Sánchez Almeida, J., Cabello, I., & Domingo, V. 2008, ApJ, 687, L131 [NASA ADS] [CrossRef] [Google Scholar]
 Born, M., & Wolf, E. 1980, Principles of Optics: Electromagnetic Theory of Propagation of Light (Oxford: Pergamon Press) [Google Scholar]
 Braat, J., Dirksen, P., & Janssen, A. J. E. M. 2002, J. Opt. Soc. Am. A, 19, 858 [NASA ADS] [CrossRef] [Google Scholar]
 Criscuoli, S., del Moro, D., Bonet, J. A., & Márquez, I. 2005, Sol. Phys., 228, 177 [NASA ADS] [CrossRef] [Google Scholar]
 Denker, C., Yang, G., & Wang, H. 2001, Sol. Phys., 202, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Gonsalves, R. A., & Chidlaw, R. 1979, in SPIE Conf. 207, ed. A. G. Tescher, 32 [Google Scholar]
 Janssen, A. J. E. M. 2002, J. Opt. Soc. Am. A, 19, 849 [NASA ADS] [CrossRef] [Google Scholar]
 Löfdahl, M. G., & Scharmer, G. B. 1994a, A&AS, 107, 243 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Löfdahl, M. G., & Scharmer, G. B. 1994b, A&AS, 107, 243 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Löfdahl, M. G., Berger, T. E., Shine, R. S., & Title, A. M. 1998, ApJ, 495, 965 [NASA ADS] [CrossRef] [Google Scholar]
 Mathar, R. J. 2009, Serbian Astron. J., 179, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Nijboer, B. R. A. 1942, Ph.D. Thesis, University of Groningen [Google Scholar]
 Noll, R. J. 1976, J. Opt. Soc. Am., 66, 207 [NASA ADS] [CrossRef] [Google Scholar]
 Paxman, R. G., Schulz, T. J., & Fienup, J. R. 1992, J. Opt. Soc. Am. A, 9, 1072 [NASA ADS] [CrossRef] [Google Scholar]
 Rayrole, J., & Mein, P. 1993, in The Magnetic and Velocity Fields of Solar Active Regions, IAU Colloq. 141, ed. H. Zirin, G. Ai, & H. Wang, ASP Conf. Ser., 46, 170 [Google Scholar]
 Rimmele, T. R. 2000, in SPIE Conf. Ser. 4007, ed. P. L. Wizinowich, 218 [Google Scholar]
 Roddier, N. 1990, Opt. Engin., 29, 1174 [NASA ADS] [CrossRef] [Google Scholar]
 Scharmer, G. B., Shand, M., Lofdahl, M. G., Dettori, P. M., & Wei, W. 2000, in SPIE Conf. Ser. 4007, ed. P. L. Wizinowich, 239 [Google Scholar]
 Scharmer, G. B., Gudiksen, B. V., Kiselman, D., Löfdahl, M. G., & Rouppe van der Voort, L. H. M. 2002, Nature, 420, 151 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 van Noort, M., Rouppe van der Voort, L., & Löfdahl, M. G. 2005, Sol. Phys., 228, 191 [NASA ADS] [CrossRef] [Google Scholar]
 von der Lühe, O. 1993, A&A, 268, 374 [NASA ADS] [Google Scholar]
 von der Lühe, O. 1994, A&A, 281, 889 [NASA ADS] [Google Scholar]
 Wang, J. Y., & Markey, J. K. 1978, J. Opt. Soc. Am., 68, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Wöger, F., & von der Lühe, O. 2007, Appl. Opt., 46, 8015 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Wöger, F., von der Lühe, O., & Reardon, K. 2008, A&A, 488, 375 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Footnotes
 ...Mathar (2009)^{}
 See also http://arxiv.org/abs/0809.2368
 ...)^{}
 Some advantage is gained by using instead the KarhunenLoève expansion. The reason is that these functions are empirically built to explain the largest amount of turbulence variance with the fewer number of functions.
All Figures
Figure 1: Example of unique basis functions for representing the PSF, including terms up to the Noll index j=7. These functions are calculated for a telescope of 100 cm diameter at a wavelength of 5250 Å, a typical situation in presentday telescopes. The spatial dimensions are in units of pixels, that we choose to be of size 0.055 , also typical of present observing conditions. The upper panel shows the results for the focused image while the lower panel shows the results for the defocused image, where the defocus is identical to that shown in Sect. 5. 

Open with DEXTER  
In the text 
Figure 2: Histograms of the distribution of the wavefront expansion coefficients for an atmosphere with Kolmogorovtype turbulence (upper row). The real and imaginary part of the corresponding expansion coefficients of the field distribution are shown in the middle and lower rows, respectively. The simulations are carried out for a telescope of D=90 cm and a Fried parameter r_{0}=7 cm. 

Open with DEXTER  
In the text 
Figure 3: Reconstruction example at 395 nm. The lower image shows the infocus image while the right image presents the reconstructed image using only one pair of infocus and defocused images. The contrast of the original image is 9.6% while the reconstructed one presents a contrast of 11.9%. Data was acquired with the SST telescope. 

Open with DEXTER  
In the text 
Figure 4: Example of the phasediversity reconstruction at 850 nm. We represent the focused image (left) and the reconstructed subimage indicated with a box (right). Images were acquired with THEMIS. 

Open with DEXTER  
In the text 
Copyright ESO 2010