3.1

- Joint estimation

The general approach for jointly estimating the object and the aberrations, which is sometimes called the joint Maximum A Posteriori (JMAP) approach, is defined by

where f(i₁, i₂, o, a; θ) is the joint probability density function of the data (i₁,i₂), of the object o and of the aberrations a. It also depends on the set of hyperparameters θ (θ stands for the regularization parameters). f(i₁|o, a; θ) and f(i₂| o, a; θ) denote the likelihood of the data i₁ and i₂, f(o; θ) and f(a; θ) are the a priori probability density function of o and a. The noise in the images is the sum of the photon noise (Poisson distributed random variable) and the Gaussian detector noise. In the case of Earth scenes, the mean signal is high enough so that the photon noise can be approximed as Gaussian. Besides, we can assume that the global noise is a stationary white Gaussian noise with a variance σ² (the same for the two images). We choose a Gaussian prior probability distribution for the object with a mean of o_m and a covariance matrix R_o also a Gaussian one for the aberrations with a null mean and a covariance matrix R_a. Hence we have

where |x| denotes the determinant of x, N² is the number of pixels in the image and k the number of aberrations. The criterion to minimize, L_JMAP(o, a) = − ln(f(i₁, i₂, o, a; θ)), can be written in the Fourier domain. Furthermore its derivative with respect to the object gives[10] a closed-form expression for the object ô(a) that minimizes the criterion for a given a.

where tilde denotes the two-dimensional DFT and S_o is the power spectral density of the object. Substituting ô(a) into the criterion yields a new criterion that does not depend explicitly on the object:

The a priori information required on the object consists of the power spectral density of the object. We choose the following model

where E stands for the mathematical expectation. This heuristic model and similar ones have been quite widely used[11]. Notice that this model requires the choice of three hyperparamaters θ_o = (k, v_o, p).

We will note θ_a the hyperparameters due to the aberrations prior, and θ_n those due to the noise statistics.

3.2

- Marginal estimation

3.2.1

- Criterion

When phase diversity is used as a wavefront sensor, the aberrations only need to be estimated, the observed object is a nuisance parameter. We propose to estimate the sole phase or the phase and the hyperparameters by integrating the object out of the problem. This novel estimator is obtained by integrating the probability density function used for the joint estimation and it is based on the Maximum A Posteriori estimator for a, which is defined by

We denote I = (i₁ i₂)^T the vector which concatenates the data. I depends linearly on variables (o and n) that are Gaussian, so it follows a Gaussian statistics. Maximizing f(i₁, i₂, a; θ) = f(I, a; θ) is thus equivalent to minimizing the following criterion:

where m_I = (H₁o_m H₂o_m)^T is the vector which contains the mean images, R_I is the covariance matrix of I and |R_I| its determinant. R_I has the following expression (1 is the identity matrix):

3.2.2

- Relationship between the joint estimator and the marginal estimator

It can be shown[12] that the two criterions L_MAP and L′_JMAP are related by the following relationship:

Thus, the difference between the marginal estimator and the joint estimator consists of a single additional term, which in the Fourier domain is given by:

Although the two estimators differ only by one term, we shall see in section 4 that their properties differ considerably. Before that, let us see how the hyperparameters can be estimated.

3.3

- Estimation of the hyperparameters

The marginal estimator as the joint one calls for the choice of the values of the hyperparameters θ = (θ_n, θ_o, θ_a). The difficulty for estimating θ_n, θ_o and θ_a is not equal : the hyperparameters of the aberrations θ_a could be easily adjusted by using Kolmogorov statistics, the hyperparameters of the object θ_o must be estimated for each object because they depend on the structure of the object. For the marginal estimator, the estimation of θ_o and θ_n is not a problem because they can be easily estimated jointly with the aberrations. Indeed, they are given by

There are four hyperparameters θ_n = σ², θ_o = (k, v_o, p) to estimate but by doing the following change of variable μ = σ²/k, we can find a closed-form expression for such that = 0. The criterion to minimize becomes L_MAP(a, μ,v_o, p, θ_a). On the contrary, for the joint estimator, hyperparameters θ_o and θ_n can not be jointly estimated with o and a. Indeed, the criterion degenerates when one seeks θ_o together with o and a; in particular it is possible to find a pair (θ_oe, o_e) that does not depend on the data such that the criterion tends to minus infinity. So before minimizing L_JMAP, these hyperparameters must be chosen by the user.

4.1

- Image simulation

In this section, we shall compare the two estimators by means of simulations. The simulations have been obtained in the following way: our object is an Earth view. The phase is generated with the first 21 Zernike polynomials (their values are listed in Table 1) and the estimated phase will be expanded on the same polynomials. Note that the coefficients a_1–3 will not be estimated: the piston coefficient a₁ is a constant added to the phase and has no influence on the point-spread function and the tilts coefficients a₂, a₃ introduce a shift in the image that is of no importance for extended objects. The defocused amplitude for the second observation plane is 2π rad peak to peak. The simulated images are monochromatic and are sampled at the Shannon rate. They are obtained by the convolution (made by FFT) of the point-spread function and the object (see Figure 2) and corrupted by a stationary white Gaussian noise.

Figure 2:

(a) aberrated phase (λ/7rms),(b) true object, (c) focused image, (d) defocused image

Table 1:

Values of the coefficients used for simulations.

4.2

- Properties

One acknowledged problem with joint estimation is that it does not ensure that bias and variance asymptotically vanish, i.e. the estimate may not converge towards the true value as the data size tends to infinity[13, 14]. On the contrary, the marginal estimator, under not very restrictive conditions, has a null asymptotic bias and variance[15, 16]. Figure 3, obtained with simulations, illustrates these asymptotic behaviors. Notice that we use a gradient-conjugate method to minimize the criterion and that these curves, like all the following ones, are the average on 20 realizations of noisy monochromatic image pairs. Furthermore we have not regularized the aberrations (an explicit regularization is not necessary when the phase is expanded on few Zernike polynomials (21)). For the two estimators, we have set the value of the hyperparameters (noise + object): σ² is equal to the average number of photons per pixels and we fit the power spectral density model with the true object to obtain the hyperparameters (k, v_o and p), which will call the “true” hyperparameters. For the marginal estimator (left figure 3), the root mean square error (RMSE) of the aberration decreases when the number of data increases. On contrary, for the joint estimator (right figure), the increase of the number of data has no effect on the RMSE. An intuitive explanation of this phenomenon is that, if a bigger image is used to estimate the aberrations, the size of the object, which is jointly reconstructed, increases too, so that the ratio of the number of unknowns to the number of data does not tend towards zero. Additionally the joint criterion may present local minima, whereas according to a Central Limit Theorem, when the number of data (i.e. the image size) increases, the shape of the likelihood function tends to become more and more Gaussian so that it presents less and less local minima.

Figure 3:

RMSE of aberrations estimates as a function of noise level. Left figure is for the marginal estimator, right figure is for the joint estimator.

4.3

- Influence of the hyperparameters

The problem of the hyperparameters is not the same as we use the joint estimator or the marginal one. For the joint estimator, we have pointed out that they must be estimated by hand, we will study now the influence of the global hyperparameter (1/k) i.e. the one that quantifies the trade-off between goodness of fit to the data and the a priori, on aberrations and object estimates. Figure 4 shows the RMSE on the aberrations estimates and on object estimate as a function of the value of this hyperparameter (we varied 1/k around its true value which corresponds to 1 in figure 4). We see that the best value of this hyperparameter is not the same for the object and for the aberrations. For the estimation of the aberrations, the reconstruction is better when the object is under-regularized. The behavior of the curve corresponding to the aberrations depends on the noise level: for a high noise level (14% here), there is one optimal hyperparameter but for lower noise levels (4%), any value under 1, including a null regularization, leads to an accurate estimation of the aberrations even though the jointly estimated object is totally noisy. This observation sheds some light on the fact that, for low noise level, unregularized phase estimation has been successfully used in the literature [17, 18] and also that explains the lack of regularization used on the object[19, 2, 3, 6] (it is difficult to regularize the object without losing on the aberrations estimates). This observation also leads us to study the asymptotic behavior of the joint estimator with no regularization. Figure 6 shows the results. In this case, when the number of data increases, the RMSE on the aberrations estimates decreases. The number of data samples over the number of unknowns is the same as for the estimation with the true hyperparameters but here, the estimator behaves, surprisingly, as if the object were not being estimated. For the marginal estimator, we compare the quality of the aberrations reconstruction obtained by minimizing L_MAP with the true hyperparameters and by minimizing L_MAP(a, μ, v_o, p, θ_a) which we call unsupervised estimation, for two different image sizes. From figure 5, we see that for low noise levels, the agreement is very good. For 128 × 128 pixels, it is quite good for any noise level. For 32 × 32 pixels and high noise level, the reconstruction is seriously degraded because of the lack of information contained in the noisy data.

Figure 4:

Plots of RMSE for aberrations estimates (solid line) and object estimate (dashed line) as a function of the value of the hyperparameter 1/k for an image size of 3 × 32 pixels. Left figure is for a noise level of 14%, right figure is for 4%.

Figure 5:

Performance of marginal estimation with the true hyperparamaters (pluses +) and for an unsupervised estimation (square □) as measured by RMSE of aberrations estimates vs. noise level. The solid line indicates an image size of 128 × 128 pixels; and the dashed line indicates an image size of 32 × 32 pixels

Figure 6:

RMSE of joint aberrations estimates as a function of noise level for a null regularization. The solid line indicates an image size of 128 × 128 pixels; the dashed line indicates an image size of 64 × 64 pixels; and the dotted line indicates an image size of 32 × 32 pixels.

4.4

- Numerical comparison

We present results of the comparison between the joint estimator and the marginal estimator. In order to compare these estimators in a realistic way, we use the joint estimator with a null regularization and for the marginal approach, the hyperparameters are reconstructed jointly with the aberrations. We compared the root mean square error of the aberration estimates as a function of noise level for two image sizes (32×32 pixels and 128×128 pixels). The results for the two estimators are plotted in figure 7. We can see at least two different domains: when the noise level n < 5%, the two estimators give approximately the same results. For 5% < n < 20%, the MAP estimation is better. We have noticed that the results, for very high noise level n > 15% and an image size of 32 × 32 pixels, depends on the value of the aberrations coefficients i.e. the joint estimator could give better results than the marginal estimator. Anyway, in practice, this noise level is never reached. In terms of computational time, these two methods are quite equivalent.

Figure 7:

RMSE of unsupervised marginal estimator (square □) and the joint one with null regularization (pluses +) aberrations estimates as a function of noise level. The solid line indicates an image size of 128 × 128 pixels; and the dashed line indicates an image size of 32 × 32 pixels