Modal-based phase retrieval using Gaussian radial basis functions

P. J. Piscaer; A. Gupta; O. Soloviev; M. Verhaegen

doi:10.1364/JOSAA.35.001233

1. INTRODUCTION

In the phase retrieval (PR) problem, the phase of a complex-valued function is recovered from measurements of the magnitude of its Fourier transform. This inverse problem has many different applications in optical imaging; see Ref. [1] for a contemporary overview. Algorithmic PR-based optical wavefront reconstruction offers an attractive means of estimating the complex generalized pupil function (GPF) from a set of measurements of the point-spread functions (PSF) in adaptive optics due to its experimental simplicity [1,2]. Moreover, other additional optical components, such as beam splitters and wavefront sensors, are not necessary, avoiding problems related to non-common-path errors and a loss in observed light intensity.

PR algorithms can be divided into two subcategories. The classical and still most implemented class of algorithms is the alternating projection (AP) methods, pioneered by Gerchberg and Saxton [3] and Fienup [4]. Recently, optimization-based algorithms representing PR as a matrix completion problem have been developed [5–7]. The solution requires solving a matrix rank minimization problem, which is non-deterministic polynomial-time (NP)-hard. A convex relaxation was proposed using the trace norm as a convex surrogate to the rank operator, approximating the problem with a semi-definite program. More recently, another optimization-based approach to solve the PR problem was presented in Ref. [8]. This algorithm is shown to be superior in terms of computational complexity, making it more suitable for larger-scale PR problems. Both classes of algorithms use multiple images at different defocus planes in order to resolve non-uniqueness issues. This technique improves the stability by incorporating extra information in the intensity measurements [2] and is one possible implementation of the more general concept of structured illumination (see, e.g., [5]). A superiority in terms of convergence toward a unique solution and stability for noisy measurements makes the convex optimization-based approaches an interesting alternative to the AP methods [5,9]. Moreover, these algorithms solve the PR problem by explicitly minimizing a cost function, making it easier to introduce structures such as sparsity.

The standard approach in both classes of algorithms is to aim for the recovery of the complete GPF in a pixel basis, such that the measured PSF is the magnitude of the two-dimensional (2D) Fourier transform of the signal to be recovered. Exploiting the computational efficiency of the fast Fourier transform (FFT), AP methods can be implemented very efficiently. However, the large number of variables corresponding to the pixel-wise representation are problematic for the optimization-based algorithms, making them only suitable for small-scale applications. In this paper, we reduce the number of variables by modeling the GPF as a linear combination of modes as an alternative to the zonal pixel-by-pixel model. This approach was shown to be promising for the adaptive optics application [10], allowing the use of optimization-based PR on a conventional computer. A trade-off between approximation accuracy of the modal basis and computational effort defines the required number of modes to be used.

The complex-valued Zernike polynomials introduced as a consequence of the extended Nijboer–Zernike (ENZ) theory [11,12] formed the chosen basis in Ref. [10]. The major limitation of the global ENZ polynomial representation is that each term extends its influence over a circular pupil, making it inflexible with respect to the pupil geometry. The Zernike theory can be adapted to other pupil geometries. However, this requires a complex reformulation of the basis for every different pupil. Moreover, they are subject to Runge’s phenomenon, leading to oscillations on the edges of the domain. Recently, GPF approximation based on Gaussian radial basis functions (GRBFs) was used for semi-analytic evaluation of the diffraction integral as an alternative to ENZ polynomials. An improvement in terms of complexity, accuracy, and execution time was achieved [13]. An important feature of the GRBF is the almost local character of each function. Since the width and location of the GRBF are free to choose for each basis function, it offers a more intuitive basis to represent the GPF that is easier to implement for any arbitrary aperture geometry. This creates an increased flexibility to the geometry of the pupil function and the possibility to model local details and sharp features compared to ENZ polynomials.

This paper is concerned with the application of GRBFs as a modal decomposition of the GPF as an alternative for ENZ polynomials. The choice of several hyper-parameters that define the shape and placement of the GRBF are investigated. The relation between the numerical conditioning and the accuracy of the solution to PR problems is important in practical implementations [2,14]. When the standard representation is ill-conditioned, algorithms such as RBF-QR [15] have been proposed to transform the GRBFs into a well-conditioned basis. Therefore, this aspect of conditioning is not evaluated in this paper. To test the performance of the new method, PR simulations are performed. Two types of aberrations are considered. First, aberration data is generated from a Zernike polynomial basis with its coefficients sampled from an assumed distribution based on empirically determined correction capabilities of a deformable mirror (DM) [10]. Second, aberrations corresponding to the correction capabilities of a higher-order DM are derived experimentally to create a phase disturbance with higher spatial frequencies. The PSF is simulated at multiple planes along the optical axis, introducing phase diversity in terms of defocus. The two different classes of basis functions, GRBFs and ENZ polynomials, are compared in terms of their theoretical fitting accuracy and performance in modal-based PR algorithms.

The structure of the paper is as follows. The formulation of the modal-based PR problem as an optimization problem is presented in Section 2. An overview of the different basis functions used to approximate the GPF is also contained in this section. The aberration data generation and simulation experiment design are discussed in Section 3. The theoretical GPF fitting accuracy for the GRBFs, including the tuning of the hyper-parameters, is explained in Section 4. The simulation results for aberration retrieval for a number of different cases are presented in Section 5. Finally, the conclusions are drawn in Section 6.

2. MODAL-BASED PHASE RETRIEVAL

A mathematical formulation of the PR problem is briefly presented in this section. The effects of aberrations on an optical system can be modeled using the GPF. The GPF is a complex function [2]:

P (ρ, θ) = A (ρ, θ) \exp (i ϕ (ρ, θ)),

where

A (\cdot)

and

ϕ (\cdot)

are real-valued functions that denote the amplitude apodization function and phase aberration, respectively, and

(ρ, θ)

are the normalized polar coordinates on the exit pupil plane. Under the assumption of purely phase aberrated systems with circular exit pupils,

A (ρ, θ)

is modeled as a characteristic function of unity values inside the pupil and zero outside. The field in the focal plane is related to that in the exit pupil by the following integral:

U (r, ϕ, f) = \frac{1}{π} \int_{0}^{1} \int_{0}^{2 π} \exp (i f ρ^{2}) P (ρ, θ) \times \exp (i 2 π r ρ \cos (θ - ϕ)) ρ d ρ d θ,

where

(r, ϕ)

are the polar coordinates in the focal region normalized with respect to the axial diffraction unit (

λ / NA

), NA being the image-side numerical aperture of the optical system. The defocus parameter

f

is used to deliberately introduce a defocus aberration to the GPF and is necessary for the convergence of optimization-based PR algorithms [5].

U (r, ϕ, f)

is the complex-valued PSF corresponding to the GPF. Only the intensity image of

U (r, ϕ, f)

, called the PSF, is observed by the camera:

y (r, ϕ, f) = {| U (r, ϕ, f) |}^{2} .

In this paper, we will define the PR problem as recovering the phase aberration

ϕ (ρ, θ)

from multiple focal-plane intensity measurements

y (r, ϕ, f)

with different introduced defocus

f

. Often, we will adopt a sampled representation in which both the GPF and (complex) PSF are sampled on an equally spaced square grid, denoted by

P \in C^{N_{p} \times N_{p}}

and

U_{i} \in C^{N_{u} \times N_{u}}

, respectively, where the subscript

i

indicates the image at focal position

f = f_{i}

. The intensity measurements of

U_{i}

are vectorized into

y_{i} \in R^{N_{u}^{2}}

. The number of diversity images will be denoted by

N_{f}

.

A more concrete formulation of the PR problem requires a convenient and systematic parameterization of the GPF. The most flexible parameterization is a pixilation of the pupil, as it can be used with any pupil geometry. However, the pixel basis requires a large number of parameters to be identified using PR. The large number of variables is problematic for the optimization-based algorithms, since they cannot exploit the computational efficiency of the FFT, as is done by AP algorithms. Parameterizations based on approximating the GPF as a linear superposition of a small number of basis functions reduces the size of the problem dramatically, since it requires estimation of just a complex scalar coefficient for each of the basis functions. The introduction of this modal decomposition allows the use of the computationally demanding optimization-based algorithms on a conventional computer [10]. Next, two different modal representations are presented.

A. Extended Nijboer–Zernike Polynomials

The representation of phase aberration $ϕ$ in terms of Zernike polynomials was generalized to represent the GPF under the ENZ theory [11]. The GPF is approximated as a truncated series of ENZ polynomials [12],

{\hat{P}}_{E} (ρ, θ) = \sum_{n, m} β_{n}^{m} N_{n}^{m} (ρ, θ) .

n

and

m

denote, respectively, the radial order and azimuthal frequency of the ENZ polynomial

N_{n}^{m} (ρ, θ)

(see Appendix A). The polynomials are ordered according to their radial order, such that the coefficients can be collected into a single vector

β \in C^{N_{β}}

, where

N_{β} = (n_{M} + 1) (n_{M} + 2) / 2

,

n_{M}

being the maximum radial order considered.

B. Gaussian Radial Basis Functions

Alternatively, the pupil function can be approximated by a linear combination of GRBFs [13]. The complex GPF is approximated by a real-valued, radially symmetric GRBF,

{\hat{P}}_{R} (ρ, θ) = A (ρ, θ) \sum_{k = 1}^{N_{γ}} γ_{k} Ψ_{k} (ρ, θ),

Ψ_{k} (ρ, θ) = e^{- λ_{k} (ρ^{2} + ϱ_{k}^{2} - 2 ρ ϱ_{k} \cos (θ - ϑ_{k}))},

where

γ_{k} \in C

,

(ϱ_{k}, ϑ_{k})

are the polar coordinates of the GRBF nodes on a polar grid and

A (ρ, θ)

is the same as in Eq. (1). Also,

λ_{k} > 0

is the shape parameter inversely proportional to the width of the GRBF. Numerical conditioning of the basis is an important issue in RBF approximation. It should be noted that, in contrast to ENZ polynomials, the GRBFs are not orthogonal. Severe ill-conditioning can occur in the flat basis function limit

(λ_{k} \to 0)

. In the literature, methods have been proposed that yield a well-conditioned basis, in which the basis functions are different, but together span the same space as the original RBF set (see, e.g., Ref. [15]). However, for the applications in this paper,

λ_{k}

was chosen to be large enough so that no problems occurred due to ill-conditioning.

C. Phase Retrieval as an Optimization Problem

Both modal representations in Eqs. (4) and (5) can be expressed in the following form:

\hat{P} (ρ, θ) = \sum_{k = 1}^{N_{α}} α_{k} B_{k} (ρ, θ)

with coefficients

α_{k} \in C

and the basis functions

B_{k} (ρ, θ)

represent either the ENZ polynomials or GRBFs. By sampling the GPF on a regular

N_{p} \times N_{p}

grid, such that each pixel corresponds to a location

(ρ_{i}, θ_{i}), i = 1, \dots, N_{p}^{2}

, the basis functions and the estimated GPF can both be represented as a matrix:

\hat{P} = \sum_{k = 1}^{N_{α}} α_{k} B_{k}

with

B_{k} \in C^{N_{p} \times N_{p}}

and

\hat{P} \in C^{N_{p} \times N_{p}}

. After vectorizing this representation, the modal decomposition becomes a single matrix-vector multiplication:

\hat{p} = B α,

such that

\hat{p} \in C^{N_{p}^{2}}

,

α \in C^{N_{α}}

, and the

k

th column of

B \in C^{N_{p}^{2} \times N_{α}}

is a vectorized representation of

B_{k}

.

Due to the linearity property of the diffraction integral in Eq. (2), the predicted PSF is a linear combination of transformed basis functions weighted by the same coefficients as the GPF. This transformation is performed using the 2D discrete Fourier transform (DFT) denoted by $F {\cdot}$ . A sampled representation of $U$ is represented on an $N_{u} \times N_{u}$ grid, where $N_{u} = D N_{p}$ , $D \geq 1$ being a constant related to the diffraction limit of the optical system. The increase in dimensions $N_{u} - N_{p}$ is computationally realized by zero-padding of the GPF before taking the 2D-DFT. For an image at a position along the optical axis corresponding to a defocus parameter of $f_{d}$ , the estimated complex image ${\hat{U}}_{d}$ can be written as

{\hat{U}}_{d} = \sum_{k = 1}^{N_{α}} α_{k} C_{d, k},

where

C_{d, k} = F {B_{k} ⊙ P_{d}} \in C^{N_{u} \times N_{u}}

,

⊙

represents the element-wise product, and

P_{d}

is the defocus function

\exp (i f_{d} ρ^{2})

sampled on the same

N_{p} \times N_{p}

grid as the basis functions extended with the correct zero-padding. Also, the modal decomposition of the complex image

U_{d}

can be represented as a matrix-vector multiplication, i.e., we define the vectorization of

{\hat{U}}_{d}

as

{\hat{u}}_{d} \in C^{N_{u}^{2}}

, such that

{\hat{u}}_{d} = C_{d} α

with the

k

th column of

C_{d} \in C^{N_{u}^{2} \times N_{α}}

containing the vectorized representation of

C_{d, k}

. Finally, the estimated intensity measurement (PSF) is given by

{\hat{y}}_{d} = {| C_{d} α |}^{2} .

The PR problem is formulated in this paper as the minimization of the error between the measurements $y_{d}$ and the estimated PSF ${\hat{y}}_{d}$ , leading to the following optimization problem:

\min_{α \in C^{N_{α}}} \sum_{d = 1}^{N_{f}} ‖ y_{d} - {| C_{d} α |}^{2} ‖,

where

‖ \cdot ‖

denotes a vector norm of interest. This problem will be solved using an efficient optimization-based PR algorithm called convex optimization-based phase retrieval (COPR) [8]. In principle, other optimization-based algorithms such as PhaseLift can be used. However, they are too computationally demanding for even medium sized problems. Since the goal of this paper is to show the advantage of GRBFs for high-order aberrations, which require a larger number of basis functions, COPR is chosen because of its superior speed.

3. SIMULATION DESIGN

This section discusses the simulations that will be performed to analyze the advantage of using GRBFs over ENZ polynomials. By considering their advantages and disadvantages, as outlined in the introduction, we expect to show an improvement of GRBFs over ENZ polynomials for several cases. First of all, GRBFs are per definition more flexible to different pupil geometries, since ENZ polynomials are defined over the unit disk only. Moreover, it is expected that the GRBF representation is more suitable to fit higher spatial frequencies in the GPF. To validate this, we will test the PR problem on both low-order and high-order aberrations. The generation of both phase types is discussed below. Finally, the implementation details for the performed simulation experiments are presented.

A. Generation of Low-Order Aberrations

For the generation of aberrations containing only lower spatial frequencies, the phase is represented in terms of Zernike polynomials:

ϕ (ρ, θ) = \sum_{n, m} ζ_{n}^{m} Z_{n}^{m} (ρ, θ),

where the Zernike polynomials

Z_{n}^{m}

are defined in Appendix A, and the coefficients

ζ_{n}^{m}

are drawn from an assumed distribution. The distribution is based on the experimentally derived correction capabilities of a low-order membrane DM [10], which followed approximately an exponential decrease in the values for

ζ_{n}^{m}

for increasing radial orders. The total number of Zernike terms considered is denoted by

N_{z} = (n_{z} + 1) (n_{z} + 2) / 2

, where

n_{z}

is the maximum radial order considered. In the following,

N_{z} = 66

, i.e., we consider Zernike polynomials up to the tenth radial order. Moreover,

ζ_{0}^{0}, ζ_{1}^{- 1}, ζ_{1}^{1} = 0

, i.e., the piston and tip–tilt modes are not included. The remaining 63 coefficients are vectorized after ordering them by increasing the radial order (Noll’s sequential index) into

ζ

. The standard deviation of the normal distribution generating the

k

th index

ζ_{k}

is computed as

c_{l} \exp (k)

, where

c_{l}

is a scaling factor to control the amplitude of the aberration. An example of such a generated wavefront is shown in Fig. 1(a).

Fig. 1. Example of a low-order and high-order aberration generated as described in Sections 3.A and 3.B, respectively. The amplitude of the phase can be scaled to any desired value.

Download Full Size | PDF

B. Generation of High-Order Aberrations

For the generation of aberrations with higher spatial frequencies, a similar experiment as in Ref. [10] is repeated for a high-order DM, having $N_{m} = 952$ actuators (Boston Micromachines KiloDM [16]). The control signal to each of the actuators $u_{i}$ is collected into the vector $u \in R^{N_{m}}$ . Because of the large number of actuators, an accurate Zernike representation of the wavefront would require too many parameters. Therefore, the DM influence matrix $H \in R^{N_{s} \times N_{m}}$ is identified such that it satisfies the relation $s = H u$ , where $s \in R^{N_{s}}$ contains all the local slopes in the wavefront measured with a high-resolution Shack–Hartmann wavefront sensor ( $N_{s} = 42632$ ). For the identification of matrix $H$ , a set of 2000 random vectors ${\tilde{u}}_{1}, \dots, {\tilde{u}}_{2000}$ are drawn from a normal distribution, and corresponding measurements ${\tilde{s}}_{1}, \dots, {\tilde{s}}_{2000}$ are collected. The matrix $H$ is identified, enforcing a sparsity pattern in the matrix, such that the number of non-zero elements is 14.01% of the total number of elements. The accuracy of the identification is very sensitive to the measurement noise of the sensor and optical misalignments in the experimental setup. After a careful calibration, the identification of $H$ resulted into an average variance accounted for (VAF) [17], defined as

VAF = (1 - \frac{\sum_{i = 1}^{2000} {‖ {\tilde{s}}_{i} - H {\tilde{u}}_{i} ‖}_{2}^{2}}{\sum_{i = 1}^{2000} {‖ {\tilde{s}}_{i} ‖}_{2}^{2}}) \cdot 100 %

of 84.31%. A set of sensor measurements

s_{i}

is simulated using random input vectors

u_{i}

drawn from a normal distribution, such that

s_{i} = c_{h} H u_{i}

,

c_{h}

being a scaling parameter to control the phase amplitude. The corresponding wavefront

ϕ_{i}

, sampled on a regular

N_{p} \times N_{p}

grid, is reconstructed from

s_{i}

using the method described in Ref. [18]. The resolution

N_{p}

is limited by the resolution of the Shack–Hartmann sensor at

N_{p} \leq 146

. In this paper,

N_{p} = 64

is used to simulate the aberration unless mentioned otherwise. A square of

N_{p} \times N_{p}

is cut from the center of the original

146 \times 146

image. An example of a high-order wavefront aberration is shown in Fig. 1(b).

C. Approximating the Generalized Pupil Function

To quantify the theoretical accuracy for each basis to fit the GPF, the least-squares error approximating a given GPF is considered. The aberration $ϕ$ is generated as discussed in Sections 3.A and 3.B on an $N_{p} \times N_{p}$ grid. Similar to Eq. (9), a sampled representation of the true GPF can be defined as a vector $p \in C^{N_{p}^{2}}$ . By solving the following complex least-squares problem,

{\hat{α}}_{L S} = \arg \min_{α \in C^{N_{α}}} ‖ p - B α ‖_{2}^{2},

a theoretical best estimate of the GPF using the defined basis can be computed as

{\hat{p}}_{L S} = B {\hat{α}}_{L S}

. The normalized real-valued root mean square (RMS) of the complex-valued approximation error

p - {\hat{p}}_{L S}

, denoted by

ϵ_{p, L S}

, defines the measure of its ability to fit the GPF, i.e.,

ϵ_{p, L S} = \frac{{‖ p - {\hat{p}}_{L S} ‖}_{2}}{{‖ p ‖}_{2}} .

Since we are mainly interested in finding the correct phase

ϕ

, we define the vectorized phase as

ϕ \in R^{N_{p}^{2}}

and the phase corresponding to the estimated GPF as

{\hat{ϕ}}_{L S} \in R^{N_{p}^{2}}

, such that the normalized RMS phase error is

ϵ_{ϕ, L S} = \frac{{‖ ϕ - {\hat{ϕ}}_{L S} ‖}_{2}}{{‖ ϕ ‖}_{2}} .

Because of the dependence on the number of GRBFs

N_{α}

and the stochastic nature of the aberration, a Monte Carlo experiment is performed for various values of

N_{α}

by using 100 draws of the simulated GPF.

Besides using this experiment to show the theoretical accuracy to fit the GPF, the method in this paragraph is used to tune the GRBF hyper-parameters in the next section. In principle, the RMS error of the PSF fit obtained by solving the PR problem can also be used. However, since this is much more computationally demanding than computing the least-squares fit of the GPF directly, it would be too time consuming and is therefore not considered in this paper.

D. Phase Retrieval Experiment

The PR problem defined in Eq. (13) is solved using the COPR algorithm [8]. First, the PSF data is simulated by taking the 2D-FFT of the GPF corresponding to the simulated aberration $ϕ$ as discussed in Section 2 for $D = 2$ . Four diversity images at $f_{1} = 0$ , $f_{2} = - 1$ , $f_{3} = 2$ , and $f_{4} = - 3 rad$ are computed and vectorized into $y_{i}, i = 1, 2, 3, 4$ . A Monte Carlo experiment solving the PR problem for 25 different aberration realizations is performed for a certain combination of $N_{α}$ and aberration type (low-/high-order).

From the solution of the PR problem, we obtain a solution $\hat{α}$ . An estimate of the PSF and GPF is obtained via ${\hat{y}}_{i} = C_{i} \hat{α}$ and $\hat{p} = B \hat{α}$ , respectively. The normalized RMS error (RMSE) of the PSF and phase will be used as a scalar quantity to express the accuracy. One major difficulty in comparing the phase is caused by the fact that the PR solution provides an estimate up to a constant offset (piston), and each pixel is wrapped on the range $[- π, π]$ . This is solved by extracting the phase of the estimated GPF, unwrapping it using a 2D phase unwrapping algorithm [19], and removing the piston. Similar to Eq. (18), the real-valued normalized RMSE of the phase is defined as

ϵ_{ϕ} = \frac{{‖ ϕ - \hat{ϕ} ‖}_{2}}{{‖ ϕ ‖}_{2}} .

The cost function in Eq. (13) minimizes the PSF fitting error rather than the GPF. A measure representing the fit of the PSF is defined in a similar way:

ϵ_{y} = \frac{\sum_{i = 1}^{N_{f}} {‖ y_{i} - {\hat{y}}_{i} ‖}_{2}}{\sum_{i = 1}^{N_{f}} {‖ y_{i} ‖}_{2}} .

The results of this experiment are discussed in Section 5.

4. FITTING ACCURACY OF GRBFs

The increased flexibility of the GRBF follows directly from the introduced freedom in terms of its hyper-parameters. Each single GRBF has two parameters: its center location pair $(ϱ_{k}, θ_{k})$ and shape parameter $λ_{k}$ . Together they define the shape and location of the GRBF in the pupil plane. If they would be chosen independent from each other, they would introduce $2 N_{α}$ hyper-parameters to the PR problem. Estimating them is a highly non-linear problem usually solved by cross-validation [13]. In this section, by using their physical interpretation in the imaging application, we propose to reduce the number of parameters to one single shape parameter $λ$ and a predefined node distribution.

A. Node Distribution

Instead of choosing the center $(ϱ_{k}, ϑ_{k})$ of each basis function separately, a number of fixed configurations are considered. In this way, the centers are no longer a hyper-parameter to be determined. Since we assume to have no specific a priori knowledge of the GPF, we are looking for a general distribution that is able to fit any generic GPF realization. There are many regular and quasi-random grids that can work as a suitable node distribution. Examples are a rectangular grid with equally spaced points, a Halton-points-based grid generated from quasi-random number sequence [20], and a planar Fibonacci grid defined using a spiral represented in polar coordinates, i.e., for the $k$ th point [14],

ϱ_{k} = ϱ_{0} \sqrt{k - 1 / 2}, ϑ_{k} = 2 π k / ϕ,

where

ϱ_{0}

is an arbitrary scale factor, and

ϕ = (1 + \sqrt{5}) / 2

is the Golden ratio. These distributions are visualized in Fig. 2. All grids are defined over an area slightly larger than the unit disk (a disk with radius 1.05) to deal with the Gibbs phenomenon [13].

Fig. 2. Examples of node distributions on a 2D grid for a unit disk pupil aperture: (a) rectangular, (b) Halton, and (c) Fibonacci.

Download Full Size | PDF

The Fibonacci grid has been proved to be a competitive and robust choice when the shape parameter is optimally chosen [14]. For completeness, the simulations in this paper have been implemented for all three different grids in Fig. 2. This showed a performance that is similar for all grids, but slightly in favor of the Fibonacci configuration. Therefore, only the Fibonacci configuration is included in the results of this paper.

B. Shape Parameter

The choice of shape parameter is significant, as it affects numerical stability, accuracy of fit, and speed of convergence. The practical design of the shape parameter is data dependent, in that it depends on the variance of wavefront aberration and its spatial frequency content as well. As the data of the GPF is not available beforehand, it is desirable to find a systematic empirical approach of shape parameter selection. To reduce the number of hyper-parameters, a single constant $λ$ is assumed for each GRBF. Moreover, especially for the regular grids (rectangular and Fibonacci), a constant $λ$ is a reasonable assumption when there is no a priori knowledge of the GPF.

To determine the value of the parameter $λ$ , we compare the RMS of the approximation error for different values of $λ$ with the method discussed in Section 3.C. First of all, since the selection of the hyper-parameter is dependent on the chosen node distribution and the total number of nodes, the influence of $λ$ is shown for different values of $N_{α}$ . Moreover, as it is shown to depend on the type of aberration, the test is repeated for both the low-order and high-order aberrations discussed in Sections 3.A and 3.B. Figure 3 shows $ϵ_{p, L S}$ , as defined in Eq. (17), for various combinations of $λ$ and $N_{α}$ on the Fibonacci grid. It is clear that both for the low-order and high-order aberration, the trend is very similar. Both start with a steep increase in performance by increasing $λ$ up to a value that is dependent on $N_{α}$ . At a certain point, which we will denote by $λ^{*}$ , the steep descend turns into a more constant and smoother curve. This point, indicated by the asterisks in Fig. 3, is the same for both the low-order and high-order aberrations and seems to be depending only on $N_{α}$ . However, where the low-order aberration $ϵ_{p, L S}$ starts to increase immediately, it still slightly decreases for the high-order aberration. As can be seen in Fig. 3, the optimal value of $λ$ , $λ_{opt}$ , is in the range of 1 to 10 for the low-order and between 10 and 40 for the high-order aberrations.

Fig. 3. Influence of $λ$ on the mean value of $ϵ_{p, L S}$ in a Monte Carlo simulation. The same simulation is performed on both low- and high-order aberration data discussed in Sections 3.A and 3.B for $N_{p} = 64$ on a circular aperture with a radius of 1.

Download Full Size | PDF

In practice, the nature of the GPF is unknown, and one does not have a set of GPF with the same statistical properties as the one to be estimated. This makes the optimal value of $λ$ more difficult to find. However, $λ^{*}$ can be used as a lower bound to the selection of $λ$ , and it can be derived from knowing only the GRBF center locations. It is shown that $ϵ_{p, L S}$ is not very sensitive for choosing $λ$ larger than $λ^{*}$ . Therefore, it is a safe choice to select $λ > λ^{*}$ , depending on the desired smoothness of the reconstruction of the wavefront.

As pointed out in Section 2, the basis can become severely ill-conditioned when $λ \to 0$ . This should be kept in mind when selecting the value of $λ$ when the number of the basis function is small. In this case, it is advised to choose $λ ≫ λ^{*}$ to avoid problems while solving the PR problem. One can use the RBF-QR algorithm [15] to improve the conditioning and transform the obtained coefficients back into the original basis. As the focus of this paper is on the high-order aberrations with a relatively large amount of basis functions, this has not been implemented. The value of $λ$ during the PR simulations was increased manually when the found $λ_{opt}$ caused problems regarding ill-conditioning.

C. Comparison to ENZ Polynomials

After $λ_{opt}$ is selected for the chosen node distribution, the fitting accuracy can be compared to that of ENZ polynomials. The mean phase approximation errors $ϵ_{ϕ, L S}$ for both the low-order and high-order aberrations are summarized in Table 1. From Table 1, we can conclude that the mean values of $ϵ_{ϕ, L S}$ indicate a higher fitting accuracy of GRBF over ENZ polynomials on average. It should be mentioned that the variance over the Monte Carlo draws was found to be approximately within 10% of the mean values for high-order aberrations, but was much more significant for low-order aberrations with values of the same order of magnitude as its mean. A more comprehensive analysis discussing the importance of this large variance is included in Section 5. Moreover, the final value of $ϵ_{ϕ, L S}$ for the high-order aberration is still significant with a minimum of 0.11 for $N_{α} = 275$ . This should be kept in mind when performing the PR simulations. The performance shown here gives a theoretical minimum of the error for each specific set of basis functions. Because the PR simulation will use the RMS of the phase error [ $ϵ_{ϕ}$ in Eq. (19)] as a measure, we have only presented the theoretical minimal phase errors $ϵ_{ϕ, L S}$ is this section. The obtained errors from the PR simulation will be compared to the values in this table to validate if the algorithm has converged to the optimal solution. A similar table can be made for $ϵ_{p, L S}$ , but it is not included here since it shows an equivalent trend.

Table 1. Mean Values of $ϵ_{ϕ, L S}$ Using GRBF and ENZ Polynomials in the Monte Carlo Simulation^a

View Table

5. PHASE RETRIEVAL SIMULATION RESULTS

As discussed in Section 3, a simulation is performed to analyze the performance of the different basis functions to the PR problem. In this section, a number of cases will be considered that demonstrate the advantages of using GRBFs in modal-based PR. For each combination of $N_{α}$ and aberration type, an experiment, as described in Section 3.D, is performed to find the value of $λ$ that gives the best fit in the least-squares sense. The scaling constants $c_{l}$ and $c_{h}$ introduced in Sections 3.A and 3.B are chosen, such that the average RMS value of $ϕ$ in the Monte Carlo simulation $RMS (ϕ) \approx 0.75$ is for both low- and high-order aberrations. The performance of the method for higher values of $RMS (ϕ)$ will be briefly discussed in the following paragraph.

A. Low- and High-Order Aberrations

To show the performance of the GRBFs for aberrations containing different spatial frequencies, the PR problem is solved for both the low- and high-order aberrations. The simulation results when considering only low-order aberrations are visualized in Figs. 4(a) and 4(c). A clear decrease of the normalized RMSE is visible for both the PSF and phase, as the number of used basis functions increase. Both the ENZ polynomials and GRBFs obtain a very accurate fit, approaching the theoretical optimum of Table 1. The variance over the multiple draws in the Monte Carlo simulation is relatively large, such that both methods lead to a roughly equivalent performance.

Fig. 4. Comparison of the results of the phase retrieval simulation between lower- and higher-order aberrations for varying number of basis functions $N_{α}$ . The boxplots show the normalized RMSE $ϵ_{y}$ and $ϵ_{ϕ}$ [see Eqs. (19) and (20)] for a circular pupil with $N_{p} = 64$ , and $RMS (ϕ) \approx 0.75 rad$ . The boxes indicate the 25th and 75th percentile of the results in the Monte Carlo simulation. Lines are drawn through the medians of the data. Data outliers due to remaining phase ambiguities are discarded.

Download Full Size | PDF

The PR simulation results for the high-order aberrations are shown in Figs. 4(b) and 4(d). The GRBFs still approximate the theoretical value of Table 1 quite closely. On the other hand, the ENZ polynomials fit starts to deviate more from the theoretical minimum when more basis functions are considered. This phenomenon together with the much less significant variance than in the low-order aberration case make the differences appear much clearer. The decrease in performance from the theoretical minimum indicates that the PR algorithm is less able to converge to the optimal solution when ENZ polynomials are used.

To investigate the performance of the method for higher phase amplitudes, the PR simulations for both the low-order and high-order aberrations are repeated for a higher value of $RMS (ϕ)$ . Figure 5 shows the normalized RMS phase error $ϵ_{ϕ}$ for $RMS (ϕ) \approx 1.5 rad$ . All other conditions are kept the same as in the experiment above, such that they can be compared to Figs 4(a) and 4(b), where $RMS (ϕ) \approx 0.75$ . Both errors increase when $RMS (ϕ)$ increases. The same improving trend in the performance is visible when using more basis functions.

Fig. 5. Phase retrieval results for low-order and high-order aberrations for $RMS (ϕ) \approx 1.5 rad$ . The presentation of the results is similar to Fig. 4.

Download Full Size | PDF

From this analysis, the GRBFs appear most beneficial when approximating aberrations containing high spatial frequencies with a relatively large amount of basis functions. An example of such an approximation is shown in Fig. 6 for $N_{p} = 128$ and $N_{α} = 377$ . One thing that stands out when comparing the estimations with GRBFs and ENZ polynomials are the ringing artifacts that appear when using ENZ polynomials. Since these rings are a consequence of the PR solution and do not appear in the least-squares approximation, they cause a gap between the theoretical value of Table 1 and the obtained PR solution when using ENZ polynomials.

Fig. 6. Example of the retrieved phase for a high-order aberration, $N_{p} = 128$ , and $N_{α} = 377$ . From left to right, the first row shows the true phase aberration, the retrieved phase estimate using GRBF, and the retrieved estimate using ENZ. The second row shows the errors for GRBF ( $ϵ_{ϕ} = 0.42$ ) and ENZ ( $ϵ_{ϕ} = 0.57$ ), respectively. The figures are truncated to the color scale shown.

Download Full Size | PDF

B. Non-Circular Aperture

An important property of the GRBF is its independence on the pupil geometry. Since each basis function can be centered around an arbitrary location on the pupil, it is possible to concentrate the information on any specific shape. In contrast, ENZ polynomials are defined on only the unit circle. Although the Zernike theory can be adapted for other pupil geometries, it requires complex theoretical reformulation (see, e.g., [12]). The GRBFs provide a much more simpler implementation, in which the user can locate the basis functions. As a result, the basis can be easily adapted to any arbitrary shape without needing any complex theory.

To demonstrate the advantage of this locality, an extreme case is considered, in which an annular pupil is defined with a unity outer radius and an inner radius of 0.7. The centers of the GRBFs are located on a ring with an inner radius of 0.65 and an outer radius of 1.05 following the Fibonacci node distribution. Compared to the circular aperture, the number of non-zero pixels in the GPF has decreased. Therefore, it would be expected that when the same number of basis functions are considered, the normalized error should decrease. The PR results for this annular aperture are summarized in Fig. 7, showing that the difference in performance between GRBFs and ENZ polynomials is more significant than in the same situation with the circular aperture (recall Fig. 4). The normalized error for the GRBF has indeed decreased compared to the circular aperture. Because the ENZ polynomials have not been redefined on the new pupil, the error for the ENZ basis has increased significantly. This demonstrates the advantage that GRBFs have on non-circular apertures. An explanation for this decrease in performance when using ENZ polynomials becomes apparent when looking at an example of the retrieved phase in Fig. 8. Note how the estimate using GRBFs is more detailed, but the estimate corresponding to the ENZ polynomial basis shows that the oscillations around the edges have become more significant due to the thin aperture shape.

Fig. 7. Comparison of the results of the phase retrieval simulation between lower- and higher-order aberrations for varying number of basis functions $N_{α}$ . The boxplots show the normalized RMSE $ϵ_{y}$ and $ϵ_{ϕ}$ [see Eqs. (19) and (20)] for an annular pupil with $N_{p} = 64$ , and $RMS (ϕ) \approx 0.75 rad$ . The boxes indicate the 25th and 75th percentile of the results in the Monte Carlo simulation. Lines are drawn through the medians of the data. Data outliers due to remaining phase ambiguities are discarded.

Download Full Size | PDF

Fig. 8. Example of the retrieved phase for a high-order aberration, $N_{p} = 128$ and $N_{α} = 377$ . From left to right, the first row shows the true phase aberration, the retrieved phase estimate using GRBF, and the retrieved estimate using ENZ. The second row shows the errors for GRBF ( $ϵ_{ϕ} = 0.27$ ) and ENZ ( $ϵ_{ϕ} = 0.68$ ), respectively. The figures are truncated to the color scale shown.

Download Full Size | PDF

C. Gaussian Measurement Noise

Finally, zero-mean white Gaussian noise is added to the simulation measurements $y_{i}$ to consider the robustness of the estimation with respect to noise. Since the intensity measurements are per definition positive, negative values in the simulated measurement are truncated to zero. For the conciseness of this paper, only a single case regarding a low-order aberration is considered on a circular aperture, and $N_{α} = 65$ basis functions are used. The results are shown in Fig. 9. The noise is sampled from a single zero-mean Gaussian distribution, and the signal-to-noise ratio (SNR) is defined with respect to the average power of the image $y_{i}$ . The noise affects the estimate for both the ENZ and GRBF basis similarly for higher SNR. For an SNR of 5–10 dB, the ENZ polynomials give a normalized error higher than one, implying a completely inaccurate estimate.

Fig. 9. Normalized RMS phase error $ϵ_{ϕ}$ [defined in Eq. (19)] as a function of the SNR for low-order aberrations with an average $RMS (ϕ) \approx 0.75$ , using $N_{α} = 65$ basis functions. The black dotted line shows the error without noise.

Download Full Size | PDF

6. CONCLUSION

The problem of reconstructing phase aberrations using a modal approach for optimization-based PR algorithms has been considered in this paper. The otherwise too computationally demanding optimization-based algorithms [5,8] can be implemented on a standard desktop computer in this modal-based framework [10]. In this paper, the application of GRBFs to model the GPF has been explored as an alternative to the existing ENZ polynomials [10,11]. Because of its computational efficiency, the COPR algorithm [8] is used to solve the PR problem. One important advantage of GRBFs is the increased flexibility introduced by user-defined hyper-parameters determining the location and shape of each basis function. The number of hyper-parameters is reduced to a single parameter describing the size of a single GRBF by assuming a predefined distribution of the centers. Guidelines have been proposed to find the hyper-parameter that leads to the best fit. It was shown that the obtained basis using GRBFs is better able to approximate the GPF than ENZ polynomials. Moreover, the solution to the PR problem has been considered for both the GRBF and ENZ polynomial basis. Simulations have shown that GRBFs are significantly better to approximate aberrations that contain higher spatial frequencies. The increased flexibility of GRBFs has been demonstrated by solving the PR problem for an annular pupil. Finally, the robustness to Gaussian measurement noise was also in favor of GRBFs, showing a lower noise sensitivity.

APPENDIX A: REAL AND COMPLEX-VALUED ZERNIKE POLYNOMIALS

The phase aberration $ϕ$ can be analyzed by the orthogonal set of basis functions formed by the circle polynomials $Z_{n}^{m}$ introduced by Zernike,

ϕ (ρ, θ) = \sum_{n, m} ζ_{n}^{m} Z_{n}^{m} (ρ, θ),

where indices

n \in N_{0}

and

m \in Z

, respectively, denote the radial order and the azimuthal frequency of the Zernike polynomial

Z_{n}^{m}

, such that

n - | m | > 0

and even. The polynomials are given by the product of a radial polynomial

R_{n}^{| m |} (ρ)

and a trigonometric function

Θ_{n}^{m} (θ)

with suitable normalization

c_{n}^{m}

,

Z_{n}^{m} (ρ, θ) = c_{n}^{m} R_{n}^{| m |} Θ_{n}^{m} (θ),

where

c_{n}^{m} = {\begin{matrix} \sqrt{n + 1} & m = 0 \\ \sqrt{2 (n + 1)} & m \neq 0 \end{matrix}, Θ_{n}^{m} (θ) = {\begin{matrix} \cos (m θ) & m \geq 0 \\ - \sin (m θ) & m < 0 \end{matrix}, R_{n}^{m} (ρ) = \sum_{s = 0}^{(n - m) / 2} \frac{{(- 1)}^{s} (n - s)!}{s! (\frac{n + m}{2} - s)! (\frac{n - m}{2} - s)!} ρ^{n - 2 s} .

The GPF can be analyzed using a truncated series of ENZ polynomials [12],

N_{n}^{m} (ρ, θ) = \sqrt{n + 1} R_{n}^{| m |} (ρ) \exp (i m θ) .

Again, the coefficients can be collected in a single vector using Noll’s indexing. The normalization used here is as given in Ref. [10].

Funding

Seventh Framework Programme (FP7) (2007-2013); H2020 European Research Council (ERC) (339681).

REFERENCES

1. Y. Shechtman, Y. C. Eldar, O. Cohen, H. N. Chapman, J. Miao, and M. Segev, “Phase retrieval with application to optical imaging: a contemporary overview,” IEEE Signal Process. Mag. 32(3), 87–109 (2015). [CrossRef]

2. D. R. Luke, J. V. Burke, and R. G. Lyon, “Optical wavefront reconstruction: theory and numerical methods,” SIAM Rev. 44, 169–224 (2002). [CrossRef]

3. R. W. Gerchberg, “A practical algorithm for the determination of the phase from image and diffraction plane pictures,” Optik 35, 237–246 (1972).

4. J. R. Fienup, “Reconstruction of an object from the modulus of its Fourier transform,” Opt. Lett. 3, 27–29 (1978). [CrossRef]

5. E. J. Candes, Y. C. Eldar, T. Strohmer, and V. Voroninski, “Phase retrieval via matrix completion,” SIAM Rev. 57, 225–251 (2015). [CrossRef]

6. Y. Chen and E. Candes, “Solving random quadratic systems of equations is nearly as easy as solving linear systems,” in Advances in Neural Information Processing Systems (2015), pp. 739–747.

7. I. Waldspurger, A. d’Aspremont, and S. Mallat, “Phase recovery, maxcut and complex semidefinite programming,” Math. Program. 149, 47–81 (2015). [CrossRef]

8. R. Doelman, H. T. Nguyen, and M. Verhaegen, “Solving large-scale general phase retrieval problems via a sequence of convex relaxations,” arXiv:1803.02652 (2018).

9. H. H. Bauschke, P. L. Combettes, and D. R. Luke, “Phase retrieval, error reduction algorithm, and Fienup variants: a view from convex optimization,” J. Opt. Soc. Am. A 19, 1334–1345 (2002). [CrossRef]

10. J. Antonello and M. Verhaegen, “Modal-based phase retrieval for adaptive optics,” J. Opt. Soc. Am. A 32, 1160–1170 (2015). [CrossRef]

11. A. J. Janssen, “Extended Nijboer–Zernike approach for the computation of optical point-spread functions,” J. Opt. Soc. Am. A 19, 849–857 (2002). [CrossRef]

12. S. Van Haver, “The extended Nijboer–Zernike diffraction theory and its applications,” Ph.D. thesis (Delft University of Technology, 2010).

13. A. Martinez-Finkelshtein, D. Ramos-Lopez, and D. Iskander, “Computation of 2D Fourier transforms and diffraction integrals using Gaussian radial basis functions,” Appl. Comput. Harmon. Anal. 43, 424–448 (2017). [CrossRef]

14. M. Maksimovic, “Optical design and tolerancing of freeform surfaces using anisotropic radial basis functions,” Opt. Eng. 55, 071203 (2016). [CrossRef]

15. B. Fornberg, E. Larsson, and N. Flyer, “Stable computations with Gaussian radial basis functions,” SIAM J. Sci. Comput. 33, 869–892 (2011). [CrossRef]

16. Boston Micromachines Corporation, 2018, http://www.bmc.bostonmicromachines.com/pdf/Kilo-DM.pdf.

17. M. Verhaegen and V. Verdult, Filtering and System Identification: A Least Squares Approach (Cambridge University, 2007).

18. F. Roddier and C. Roddier, “Wavefront reconstruction using iterative Fourier transforms,” Appl. Opt. 30, 1325–1327 (1991). [CrossRef]

19. M. A. Herráez, D. R. Burton, M. J. Lalor, and M. A. Gdeisat, “Fast two-dimensional phase-unwrapping algorithm based on sorting by reliability following a noncontinuous path,” Appl. Opt. 41, 7437–7444 (2002). [CrossRef]

20. L. Kocis and W. J. Whiten, “Computational investigations of low-discrepancy sequences,” ACM Trans. Math. Softw. 23, 266–294 (1997). [CrossRef]

$N_{α}$	GRBF-LO	ENZ-LO	GRBF-HO	ENZ-HO
65	$6.56 \cdot 10^{- 3}$	$9.46 \cdot 10^{- 3}$	$5.16 \cdot 10^{- 1}$	$5.74 \cdot 10^{- 1}$
90	$2.23 \cdot 10^{- 3}$	$3.60 \cdot 10^{- 3}$	$4.09 \cdot 10^{- 1}$	$4.83 \cdot 10^{- 1}$
119	$5.82 \cdot 10^{- 4}$	$8.94 \cdot 10^{- 4}$	$3.18 \cdot 10^{- 1}$	$3.98 \cdot 10^{- 1}$
152	$1.48 \cdot 10^{- 4}$	$2.90 \cdot 10^{- 4}$	$2.42 \cdot 10^{- 1}$	$3.21 \cdot 10^{- 1}$
189	$5.04 \cdot 10^{- 5}$	$1.01 \cdot 10^{- 4}$	$1.89 \cdot 10^{- 1}$	$2.58 \cdot 10^{- 1}$
230	$1.42 \cdot 10^{- 5}$	$2.98 \cdot 10^{- 5}$	$1.42 \cdot 10^{- 1}$	$2.07 \cdot 10^{- 1}$
275	$5.84 \cdot 10^{- 6}$	$7.06 \cdot 10^{- 6}$	$1.10 \cdot 10^{- 1}$	$1.65 \cdot 10^{- 1}$

Modal-based phase retrieval using Gaussian radial basis functions

Abstract

1. INTRODUCTION

2. MODAL-BASED PHASE RETRIEVAL

A. Extended Nijboer–Zernike Polynomials

B. Gaussian Radial Basis Functions

C. Phase Retrieval as an Optimization Problem

3. SIMULATION DESIGN

A. Generation of Low-Order Aberrations

B. Generation of High-Order Aberrations

C. Approximating the Generalized Pupil Function

D. Phase Retrieval Experiment

4. FITTING ACCURACY OF GRBFs

A. Node Distribution

B. Shape Parameter

C. Comparison to ENZ Polynomials

5. PHASE RETRIEVAL SIMULATION RESULTS

A. Low- and High-Order Aberrations

B. Non-Circular Aperture

C. Gaussian Measurement Noise

6. CONCLUSION

APPENDIX A: REAL AND COMPLEX-VALUED ZERNIKE POLYNOMIALS

Funding

REFERENCES

Cited By

Figures (9)

Tables (1)

Equations (25)

Journal of the Optical Society of America A