EDP Sciences
Free access
Volume 518, July-August 2010
Herschel: the first science highlights
Article Number A6
Number of page(s) 9
Section Astronomical instrumentation
DOI http://dx.doi.org/10.1051/0004-6361/201014318
Published online 19 August 2010
A&A 518, A6 (2010)

Image reconstruction with analytical point spread functions

A. Asensio Ramos1,2 - A. López Ariste3

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

Context. The image degradation produced by atmospheric turbulence and optical aberrations is usually alleviated using post-facto 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 non-linear 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 Zernike-Nijboer 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 phase-diversity 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 diffraction-limited 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 post-facto 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 real-time. 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 real-time 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 artificially-induced 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 multi-object multi-frame 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 time-consuming 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 Nijboer-Zernike 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 phase-diversity. Finally, our conclusions are presented in Sect. 6.

2 Wavefront expansion

We define $P(\rho,\theta)$ to represent the wavefront of a point source beam in a pupil plane, where $\rho$ and $\theta$ 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)

$\displaystyle %
U(r,\varphi;f)=\frac{1}{\pi}\int_0^1\int_0^{2\pi} \rho {\rm d}\...
...rm i}f\rho^2} P(\rho,\theta) {\rm e}^{2\pi {\rm i}\rho r \cos(\theta-\varphi)},$   (1)

which is, in general, a complex quantity. In this expression, we define a coordinate system $(r,\varphi;f)$ over the image space, where $(r,\varphi)$ 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 $NA\ll 1$, without loss of generality for solar instruments. The $(r,\varphi;f)$ coordinates can be related to the geometrical coordinates (X,Y,Z) of the optical system by means of the following relationships
                              x = $\displaystyle X\frac{2\pi(NA)}{\lambda}, \quad y=Y\frac{2\pi(NA)}{\lambda}, \quad f=Z\frac{\pi(NA)^2}{\lambda}$  
r = $\displaystyle \sqrt{x^2+y^2}, \quad \tan \varphi =\frac{y}{x},$ (2)

where $\lambda$ 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 $P(\rho,\theta)$. In general, the wavefront is also a complex quantity that is usually represented in polar form: the modulation in amplitude is represented by $A(\rho,\theta)$, while the phase modulation is determined by $\Phi(\rho,\theta)$. If the pupil plane is fully transmissive (or perfectly reflective), we can assume $A(\rho,\theta)=1$ so that

$\displaystyle %
P(\rho,\theta)= A(\rho,\theta) {\rm e}^{{\rm i} \Phi(\rho,\theta)}={\rm e}^{{\rm i} \Phi(\rho,\theta)}.$   (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 $\vert m\vert \leq n$ and that n-|m| is even. Their definition is

         $\displaystyle %
Z_n^m(\rho,\theta)$ = $\displaystyle R_n^m(\rho) \cos (m\theta)$  
$\displaystyle Z_n^{-m}(\rho,\theta)$ = $\displaystyle R_n^m(\rho) \sin (m\theta)$  
$\displaystyle Z_n^0(\rho,\theta)$ = $\displaystyle R_n^0(\rho),$ (4)

where $R_n^m(\rho)$ 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 $(n, \ m)$ 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
$\displaystyle \Phi(\rho,\theta) = \sum_j a_j Z_j(\rho,\theta)= \sum_{n,m} a_n^m
Z_n^m(\rho,\theta).$     (5)

In a turbulent atmosphere of Kolmogorov type, the aberrations are characterized by coefficients aj 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:
$\displaystyle %
{\rm e}^{{\rm i} \Phi(\rho,\theta)}=1+{\rm i}\Phi(\rho,\theta)-\frac{1}{2}\Phi^2(\rho,\theta)+\ldots$     (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
$\displaystyle {\rm e}^{{\rm i} \Phi(\rho,\theta)} = 1 + {\rm i}\sum_j a_j Z_j(\rho,\theta)-\frac{1}{2} \Big[\sum_j a_j
Z_j(\rho,\theta)\Big]^2+\ldots$     (7)

For small aberrations ($a_j \ll 1$), 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 microscope-type 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 Zj Zk at second order and high-order products $Z_j \ldots Z_k$ 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

$\displaystyle Z_j(\rho,\theta)Z_{j'}(\rho,\theta)=\sum_k d_k Z_k(\rho,\theta),$     (8)

where the coefficients dk 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
$\displaystyle %
\prod_{j} Z_j(\rho,\theta)=\sum_k c_k Z_k(\rho,\theta),$     (9)

where the ck coefficients can be calculated from the dk coefficients of Eq. (8). As a consequence, the expansion of Eq. (7) can be written confidently as
$\displaystyle {\rm e}^{{\rm i} \Phi(\rho,\theta)}=1+\sum_{k=2} \beta_k Z_k(\rho,\theta),$     (10)

where the $\beta_k$ 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 $\beta_k$ that can be inferred from the coefficients aj 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 $\beta $ coefficients are approximately equal to the Zernike coefficients of the expansion of the wavefront
$\displaystyle %
{\rm Im}(\beta_k) \approx \alpha_k,$   (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 $\beta $ 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 $P(\rho,\theta)$ as a linear combination of Zernike functions with complex coefficients $\beta_k$, the field distribution of Eq. (1) simplifies to a sum of integrals of the form

                           $\displaystyle U(r,\varphi;f)$ = $\displaystyle \frac{1}{\pi}\int_0^1\int_0^{2\pi} \rho {\rm d}\rho {\rm d}\theta
{\rm e}^{{\rm i}f\rho^2} {\rm e}^{2\pi {\rm i}\rho r \cos(\theta-\varphi)}$  
    $\displaystyle + \sum_j \beta_j \frac{1}{\pi}\int_0^1\int_0^{2\pi} \rho {\rm d}\...
...{\rm e}^{{\rm i}f\rho^2} Z_j {\rm e}^{2\pi {\rm i}\rho r \cos(\theta-\varphi)}.$ (12)

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 Jacobi-Anger identity (Abramowitz & Stegun 1972)
$\displaystyle %
\int_0^{2\pi} {\rm e}^{2\pi {\rm i}\rho r\cos(\theta-\varphi)}{...
...\theta} {\rm d}\theta = 2\pi {\rm i}
{\rm e}^{{\rm i}m\varphi}J_m(2\pi \rho r),$     (13)

where the $J_m(2\pi \rho r)$ 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 ${\rm e}^{{\rm i}m\theta}$ (that contains both the $\cos$ and $\sin$ functions into its real and complex parts for an unsigned m index), we end up with
$\displaystyle %
2 \int_0^1 \rho {\rm d}\rho {\rm e}^{{\rm i}f\r...
..._0(2\pi\rho r)
+2\sum_j {\rm i}^m {\rm e}^{{\rm i}m\varphi} \beta_j V_n^m(r,f).$     (14)

The function Vnm(r,f) represents the following (complex) radial integral
$\displaystyle %
V_n^m(r,f) = \int_0^1 \rho {\rm d}\rho {\rm e}^{{\rm i}f\rho^2} R_n^{\vert m\vert}(\rho)
J_m(2\pi\rho r),$     (15)

which fulfills the condition
$\displaystyle %
[V_n^m(r,f)]^\dag = V_n^m(r,-f).$     (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:
$\displaystyle %
V_n^m(r) = (-1)^{\frac{n-\vert m\vert}{2}}\frac{J_{n+1}(2\pi r)}{2\pi r}\cdot$     (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

$\displaystyle %
U(r,\varphi;f) = 2 V_0^0(r,f)+2\sum_j {\rm i}^m {\rm e}^{{\rm i}m\varphi}\beta_j V_n^m(r,f).$     (18)

We note that, when no aberrations are present, the in-focus field distribution (at f=0) is proportional to
$\displaystyle %
V_0^0(r) = \frac{J_{1}(2\pi r)}{2\pi r}$     (19)

which is the function giving rise to the well known Airy spot for a circular non-aberrated aperture. At $f \neq 0$, the field distribution is proportional to V00(r,f), which, using the rapidly convergent series of Braat et al. (2002), can be written exactly as
$\displaystyle %
V_0^0(r,f) = {\rm e}^{{\rm i}f} \sum_{l=1}^\infty (-2{\rm i}f)^{l-1} \frac{J_{l}(2\pi r)}{(2\pi r)^l}\cdot$     (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 Vnm(r,f) defined above, and of which the (in-focus or out-of-focus) 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

$\displaystyle %     (21)

For the sake of simplicity in the presentation, we separate the analysis of the in-focus and the out-of-focus PSFs.

3.1 In-focus PSF

Using Eq. (18) with the first term included in the summation being j=0, the PSF results in the double summation
$\displaystyle %   (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 $j \neq k$. After tedious but straightforward algebra, the final expression for the point spread function is
$\displaystyle {\rm PSF}(r,\varphi)=8\sum \limits_{\tiny\begin{array}{c}j=0 \\  ...
... C_{jk} X_{n,n'}^{m,m'}(r,\varphi) + D_{jk} Y_{n,n'}^{m,m'}(r,\varphi) \right].$   (23)

The functions appearing in the linear expansion can be calculated as
$\displaystyle X_{n,n'}^{m,m'}(r,\varphi) = (-1)^{m'} C^{m,m'}(\varphi) V_n^m(r) V_{n'}^{m'}(r)$    
$\displaystyle Y_{n,n'}^{m,m'}(r,\varphi) = -(-1)^{m'} S^{m,m'}(\varphi) V_n^m(r) V_{n'}^{m'}(r).$   (24)

We note that the azimuthal and radial contributions are separate. The angular dependence is modulated by the functions
$\displaystyle C^{m,m'}(\varphi) = \cos \left[ \frac{\pi}{2} (m+m')+(m-m')\varphi \right]$      
$\displaystyle S^{m,m'}(\varphi) = \sin \left[ \frac{\pi}{2} (m+m')+(m-m')\varphi \right].$     (25)

The coefficients Cjk and Djk in the linear expansion of the PSF are related to the real and imaginary parts of products of $\beta $ coefficients
                Cjk = $\displaystyle C_{kj} = \epsilon_{jk} {\rm Re}(\beta_j \beta_k^\dag )$  
Djk = $\displaystyle -D_{kj} = \epsilon_{jk} {\rm Im}(\beta_j \beta_k^\dag ),$ (26)

where $\epsilon_{kj}=1-\delta_{kj}/2$. 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 Cjk coefficients and are azimuth-independent, since
$\displaystyle %
X_{n,n}^{m,m}(r,\varphi)$ = $\displaystyle \left[ V_n^m(r) \right]^2$  
$\displaystyle Y_{n,n}^{m,m}(r,\varphi)$ = 0. (27)

Because n-|m| 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
$\displaystyle %
\int_0^1 \int_0^{2\pi} r {\rm d}r {\rm d}\varphi {\rm PSF}(r,\varphi) = \sum_{j=0} \frac{\vert\beta_j\vert^2}{\pi(n+1)},$     (28)

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
$\displaystyle %
S = \vert\beta_0\vert^2 \left( \sum_{j=0} \frac{\vert\beta_j\vert^2}{(n+1)} \right)^{-1}\cdot$   (29)

\end{figure} 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 present-day telescopes. The spatial dimensions are in units of pixels, that we choose to be of size 0.055 $^{\prime \prime }$, 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 Out-of-focus 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 Vnm(r,f) functions in general are complex, the basis functions are slightly more elaborate
$\displaystyle %
{\rm PSF}(r,\varphi;f) = 8\sum \limits_{\tiny\begin{array}{c} j...
...jk} X_{n,n'}^{m,m'}(r,\varphi;f)
+ D_{jk} Y_{n,n'}^{m,m'}(r,\varphi;f) \right],$     (30)

                      $\displaystyle %
X_{n,n'}^{m,m'} (r,\varphi;f)$ = $\displaystyle (-1)^{m'}
\left[ C^{m,m'}(\varphi) \xi_{n,n'}^{m,m'}(r;f) +S^{m,m'}(\varphi) \psi_{n,n'}^{m,m'}(r;f) \right]$  
$\displaystyle Y_{n,n'}^{m,m'} (r,\varphi;f)$ = $\displaystyle (-1)^{m'}
\left[ C^{m,m'}(\varphi) \psi_{n,n'}^{m,m'}(r;f) -S^{m,m'}(\varphi) \xi_{n,n'}^{m,m'}(r;f) \right].$ (31)

The auxiliary functions appearing in the previous formulae are given by
$\displaystyle %    
$\displaystyle \psi_{n,n'}^{m,m'}(r;f) = {\rm Re}(V_n^m) {\rm Im}(V_{n'}^{m'}) - {\rm Im}(V_n^m) {\rm Re}(V_{n'}^{m'}).$   (32)

When the previous expressions are evaluated at the focus (f=0), because the Vnm functions are real, we find of course that $\psi_{n,n'}^{m,m'}=0$ and $\xi_{n,n'}^{m,m'}=V_n^m V_{n'}^{m'}$, thus recovering the expressions of Eq. (24). Likewise, for the terms of the summation in Eq. (30) for which j=k, we have
              $\displaystyle %
X_{n,n}^{m,m}(r,\varphi;f)$ = $\displaystyle \xi_{n,n}^{m,m}(r;f) = V_n^m(r;f) [{V_n^m(r;f)}]^\dag $  
$\displaystyle Y_{n,n}^{m,m}(r,\varphi;f)$ = 0. (33)

\end{figure} Figure 2:

Histograms of the distribution of the wavefront expansion coefficients $\alpha $ for an atmosphere with Kolmogorov-type turbulence (upper row). The real and imaginary part of the corresponding expansion coefficients $\beta $ 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 r0=7 cm.

Open with DEXTER

4 Statistical properties of the $\beta $ coefficients

In a Kolmogorov turbulent atmosphere, the coefficients $\alpha_j$ of the Zernike expansion of the wavefront follow a multivariate Gaussian distribution with a non-diagonal 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 r0. As a function of this parameter, we simulate random samples $\alpha_k$ 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 $\beta_k$ coefficients. The corresponding statistical analysis is presented in Fig. 2 for aberrations of order higher than two, excluding tip-tilt. The upper row shows that the $\alpha_k$ 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 $\beta_k$ 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 $\beta_k$ coefficients being equal to the $\alpha_k$coefficients. Higher-order corrections to the imaginary parts occur only at odd orders, so that the first correction already takes place at third-order, which often induces a relatively small change for turbulence that is not too strong. The Pearson linear correlation coefficient between $\alpha_k$ and ${\rm Im}(\beta_k)$ is larger than 0.95 for practically all orders, indicating that ${\rm Im}(\beta_k)$ 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 $\beta_k$ coefficients is modified only at second-order 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 $\beta_j$ 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 (D0) and another one with exactly the same aberrations plus a known defocus (Dk). In the standard image formation paradigm, the observed image after degradation by the atmosphere and the optical devices of the telescope can be written as
              D0 = F * P0 + N0  
Dk = F * Pk + Nk, (34)

where P0 and Pk are the PSFs of the atmosphere+telescope in the focused and defocused images, respectively. The quantities N0 and Nk are noise contributions in the image formation, each one characterized with a variance $\sigma_0^2$ and $\sigma_k^2$, 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

                                 P0 = $\displaystyle \vert{\rm FT}^{-1}\left\{A(\rho,\theta) \exp\left[{{\rm i} \Phi(\rho,\theta)}\right]\right\}\vert^2$  
Pk = $\displaystyle \vert{\rm FT}^{-1}\left\{A(\rho,\theta) \exp\left[{{\rm i} \Phi(\rho,\theta)+i\delta_d(\rho,\theta)}\right]\right\}\vert^2,$ (35)

where the defocus is introduced by adding a term of the form $\delta_d = \alpha^d_4 Z_4(\rho,\theta)$and ${\rm FT}^{-1}$ 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 $\alpha_j$ of the expansion of Eq. (5) that minimize the error metric (Paxman et al. 1992)

$\displaystyle %
L = \sum_{u,v} \frac{\vert\hat{D}_k \hat{P}_0 - \hat{D}_0 \hat{P}_k\vert^2}{(\vert\hat{P}_0\vert^2+\gamma \vert\hat{P}_k\vert^2)^{1/2}},$   (36)

where $\hat{D}_0$, $\hat{D}_k$, $\hat{P}_0$, and $\hat{P}_k$ are the Fourier transforms of D0, Dk, P0 and Pk, respectively. The summation is carried out over the full Fourier domain defined by the frequencies u and v, while $\gamma=\sigma_0^2/\sigma_k^2$. We note that the dependence of the metric on the $\alpha_j$ coefficients is highly non-linear 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 aj coefficients involves several Fourier transforms, which are of order $O(N \log N)$. 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)
$\displaystyle %
\hat{F} = \frac{\hat{D}_0 \hat{P}_0^\dag + \gamma \hat{D}_k \hat{P}_k^\dag }{\vert\hat{P}_0\vert^2+\gamma \vert\hat{P}_k\vert^2}\cdot$   (37)

The maximum-likelihood 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 $\hat H$ 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 Fourier-transforming the observed images and the basis functions X and Y, it is a matter of finding the values of the $\beta_j$ complex quantities that minimize the metric L. Another advantage is that the non-linear dependence of the metric function on the $\beta_j$ coefficients is less pronounced, thus helping the optimization routine to find the minimum.

To find the values of the coefficients $\beta_j$ 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 $\beta_j$coefficients. After some tedious but straightforward algebra, the gradient can be found to be

                               $\displaystyle %
\frac{\partial L}{\partial \beta_l^{\rm r,i}}$ = $\displaystyle 2\sum {\rm Re} \Bigg[ Q^2(\hat{D}_k \hat{P}_0 \!- \!\hat{D}_0 \ha...
...\rm r,i}} \!- \!
\hat{D}_0 \frac{\partial \hat{P}_k}{\beta_l^{\rm r,i}} \right)$  
    $\displaystyle -Q^4 \vert\hat{D}_k \hat{P}_0 - \hat{D}_0 \hat{P}_k\vert^2 {\rm R...
\frac{\partial \hat{P}_k}{\beta_l^{\rm r,i}} \hat{P}_k^\dag\right) \Bigg],$ (38)

where $Q=(\vert\hat{P}_0\vert^2+\gamma \vert\hat{P}_k\vert^2)^{-1/2}$ and
                           $\displaystyle \frac{\partial \hat{P}_0}{\beta_l^{\rm r}}$ = $\displaystyle 8 \sum_{k=0} \left( \beta_k^{\rm r} \hat{X}_{n,n'}^{m,m'}(r,\phi;0) -
\beta_k^{\rm i} \hat{Y}_{n,n'}^{m,m'}(r,\phi;0)\right)$  
$\displaystyle \frac{\partial \hat{P}_0}{\beta_l^{\rm i}}$ = $\displaystyle 8 \sum_{k=0} \left( \beta_k^{\rm r} \hat{Y}_{n,n'}^{m,m'}(r,\phi;0) +
\beta_k^{\rm i} \hat{X}_{n,n'}^{m,m'}(r,\phi;0)\right)$  
$\displaystyle \frac{\partial \hat{P}_k}{\beta_l^{\rm r}}$ = $\displaystyle 8 \sum_{k=0} \left( \beta_k^{\rm r} \hat{X}_{n,n'}^{m,m'}(r,\phi;f) -
\beta_k^{\rm i} \hat{Y}_{n,n'}^{m,m'}(r,\phi;f)\right)$  
$\displaystyle \frac{\partial \hat{P}_k}{\beta_l^{\rm i}}$ = $\displaystyle 8 \sum_{k=0} \left( \beta_k^{\rm r} \hat{Y}_{n,n'}^{m,m'}(r,\phi;f) +
\beta_k^{\rm i} \hat{X}_{n,n'}^{m,m'}(r,\phi;f)\right).$ (39)

The superindices r and i indicate the real and imaginary part of the coefficient $\beta_l$, 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 $\hat{X}$ and $\hat{Y}$ are the Fourier transforms of the functions defined in Eq. (31).

\includegraphics[width=8.25cm,clip]{14318fig3b.eps} }\end{figure} Figure 3:

Reconstruction example at 395 nm. The lower image shows the in-focus image while the right image presents the reconstructed image using only one pair of in-focus 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 field-of-view 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 $N_{\rm proc}$ processors were at hand, the computing time would be reduced by a factor $1/N_{\rm proc}$ 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 phase-diversity 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 1-m 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 $128\times128$ pixels ($\sim$ $4.5''\times4.5''$), an image of $1000\times1000$ 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 $\sim$5 min. When operating on a standard present-day 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 phase-diversity 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 phase-diversity applied to THEMIS data. The images were restored with patches of $100 \times 100$ pixels, equivalent to $7.5 ''\times 7.5''$. 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 multi-image information to the reconstruction algorithm (Criscuoli et al. 2005).

\end{figure} Figure 4:

Example of the phase-diversity reconstruction at 850 nm. We represent the focused image (left) and the reconstructed sub-image 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 Nijboer-Zernike 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 post-facto 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 wide-field 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 time-consuming 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 on-line 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 multi-image phase diversity (e.g., Paxman et al. 1992) or multi-object multi-frame 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.

Images 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 AYA2007-63881 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:

              $\displaystyle 2 \cos (m_1 \theta)\cos(m_2 \theta)$ = $\displaystyle \cos[(m_1-m_2)\theta]
+ \cos[(m_1+m_2)\theta]$  
$\displaystyle 2 \sin (m_1 \theta)\cos(m_2 \theta)$ = $\displaystyle \sin[(m_1-m_2)\theta]+\sin[(m_1+m_2)\theta]$  
$\displaystyle 2 \sin (m_1 \theta)\sin(m_2 \theta)$ = $\displaystyle \cos[(m_1-m_2)\theta]-\cos[(m_1+m_2)\theta].$ (A.1)

In other words, the result is always the sum of trigonometric functions depending on two azimuthal numbers $m_1 \pm m_2$. This is indeed the seed of two new Zernike functions with $m_3=m_1 \pm m_2$ if we are able to write the product of the radial polynomials with identical m3 azimuthal value. That is, we seek a product expansion for the radial polynomials of the form
$\displaystyle R_{n_1}^{m_1} R_{n_2}^{m_2} = \sum _{n3=m3}^{n_1+n_2}
g_{n_1,m_1,n_2,m_2,n_3,m_3} R_{n_3}^{m_3}.$   (A.2)

Every term of the series above has the same azimuthal number m3, coincident with the m3 of the trigonometric product and can therefore be combined with it to build a Zernike function Zn3m3. Mathar (2009) demonstrates that this is indeed feasible and provides the value of the gn1,m1,n2,m2,n3,m3 coefficient in the product expansion as
                                         gn1,m1,n2,m2,n3,m3 = 2(n3+1)  
  $\textstyle \times$ $\displaystyle \sum_{s_1=0}^{-a_1}\sum_{s_2=0}^{-a_2}\sum_{s_3=0}^{-a_3}
  $\textstyle \times$ $\displaystyle \prod_{j=1}^3 (-1)^{s_j}\left(\begin{array}{c}n_j-s_j \\  {s_j} \...
...y} \right)\left(\begin{array}{c} n_j-2s_j\\  {-a_j-s_j}\\  \end{array} \right),$ (3)

where a=-(n-m)/2. These coefficients play the same role as the Clebsh-Gordan recoupling coefficients in determining the product of two spherical harmonics.

Appendix B: Some properties of ${\rm X_{n,n'}^{m,m'}}$ and ${\rm Y_{n,n'}^{m,m'}}$ and analytical Strehl ratio

After interchanging the (n,m) and (n',m') pair of indices, the following expressions hold:

$\displaystyle \xi_{n',n}^{m',m}(r;f)$ = $\displaystyle \xi_{n,n'}^{m,m'}(r;f)$  
$\displaystyle \psi_{n',n}^{m',m}(r;f)$ = $\displaystyle -\psi_{n,n'}^{m,m'}(r;f)$  
$\displaystyle C^{m',m}(\varphi)$ = $\displaystyle C^{m,m'}(-\varphi)$  
$\displaystyle S^{m',m}(\varphi)$ = $\displaystyle S^{m,m'}(-\varphi).$ (B.1)

As a consequence, the functions Xn,n'm,m' and Yn,n'm,m' fulfill the following properties:
            $\displaystyle X_{n',n}^{m',m}(r,\varphi;f)$ = $\displaystyle (-1)^{m'-m} Y_{n,n'}^{m,m'}(r,-\varphi;f)$  
$\displaystyle Y_{n',n}^{m',m}(r,\varphi;f)$ = $\displaystyle -(-1)^{m'-m} X_{n,n'}^{m,m'}(r,-\varphi;f).$ (B.2)

For the Strehl ratio, it is possible to demonstrate that only the Airy part of the PSF (either in-focus or out-of-focus) contributes to the central peak of the PSF. The reason is that, evaluating Eq. (15) at r=0gives
$\displaystyle V_n^m(0,f) = \delta_{n0} \delta_{m0} \frac{1}{2f} \left[ \sin f + i(1- \cos f)\right],$   (B.3)

so that
           $\displaystyle X_{n,n'}^{m,m'}(r,\varphi;f)$ = $\displaystyle \delta_{n0} \delta_{m0} \delta_{n'0} \delta_{m'0} \frac{1-\cos f}{2f^2}$  
$\displaystyle Y_{n,n'}^{m,m'} (r,\varphi;f)$ = 0, (B.4)

which gives, for the f=0, a contribution of 1/4 at the center of the PSF.



...Mathar (2009)[*]
See also http://arxiv.org/abs/0809.2368
Some advantage is gained by using instead the Karhunen-Loè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

\end{figure} 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 present-day telescopes. The spatial dimensions are in units of pixels, that we choose to be of size 0.055 $^{\prime \prime }$, 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

\end{figure} Figure 2:

Histograms of the distribution of the wavefront expansion coefficients $\alpha $ for an atmosphere with Kolmogorov-type turbulence (upper row). The real and imaginary part of the corresponding expansion coefficients $\beta $ 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 r0=7 cm.

Open with DEXTER
In the text

\includegraphics[width=8.25cm,clip]{14318fig3b.eps} }\end{figure} Figure 3:

Reconstruction example at 395 nm. The lower image shows the in-focus image while the right image presents the reconstructed image using only one pair of in-focus 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

\end{figure} Figure 4:

Example of the phase-diversity reconstruction at 850 nm. We represent the focused image (left) and the reconstructed sub-image indicated with a box (right). Images were acquired with THEMIS.

Open with DEXTER
In the text

Copyright ESO 2010