## Abstract

The use of a spatial light modulator for implementing a digital phase-shifting (PS) point diffraction interferometer (PDI) allows tunability in fringe spacing and in achieving PS without the need for mechanically moving parts. However, a small amount of detector or scatter noise could affect the accuracy of wavefront sensing. Here, a novel method of wavefront reconstruction incorporating a virtual Hartmann-Shack (HS) wavefront sensor is proposed that allows easy tuning of several wavefront sensor parameters. The proposed method was tested and compared with a Fourier unwrapping method implemented on a digital PS PDI. The rewrapping of the Fourier reconstructed wavefronts resulted in phase maps that matched well the original wrapped phase and the performance was found to be more stable and accurate than conventional methods. Through simulation studies, the superiority of the proposed virtual HS phase unwrapping method is shown in comparison with the Fourier unwrapping method in the presence of noise. Further, combining the two methods could improve accuracy when the signal-to-noise ratio is sufficiently high.

© 2015 Optical Society of America

## 1. Introduction

The use of the spatial light modulator (SLM) as a wavefront sensor and corrector is well known. Examples include phase-shifting (PS) point diffraction interferometer (PDI) [1–5], phase retrieval [6] and shearing interferometer [7]. The SLM offers easy tunability to relevant wavefront sensor parameters and allows for a large number of applications in computational wavefrontsensing [8]. The PS PDI device is a common-path interferometer that generates a reference wavefront by diffraction from a pinhole. It is sensitive to small magnitude aberrations while being immune to external vibrations. A typical PS PDI involves the measurement of interferograms with phase shifts of 0, *π*/2 and *π*. A linear combination of these interferograms yields the wrapped phase [9]. An important ultimate step in interferometricwavefront sensors is phase unwrapping.

To achieve PS, some of the earlier works used partially transmitting diffraction gratings, polarization optics designs and liquid crystal devices [10–12]. Recently, a digital PS PDI was proposed with the help of an 8-bit transmitting SLM device from Holoeye (Model: LC 2012) [4]. The advantage of the digital PS PDI wavefront sensor is that it allows the measurement of aberrations in the optical system without the need for mechanically moving parts. The design of this wavefront sensor involves a partial transmission double square pinhole phase mask on the phase modulating SLM such that it allows PS and interference of the transmitted light. The fringe contrast can be controlled by adjusting SLM transmission, the size of the pinholes and spacing between them. The quality of the reference wavefront plays an important role and it does not compromise the quality of wavefront sensing if the pinhole size used to generate the reference wavefront is kept small enough.

However, the previous study [4] indicated that the wrapped phase obtained with a digital PS PDI is noisy and contained multiple phase singularities, making the phase unwrapping step in PS PDI crucial and necessitating iterative calculations for accurate performance. The question of successfully obtaining useful information from noisy interferograms is still open and application dependent [13]. Successfully implemented conventional phase unwrapping methods [14, 15] were found to fail consistently in the presence of noise exhibiting inconsistencies. Hence, they are ineffective in a digital PS PDI scenario. Consequently, advanced phase unwrapping methods are necessary [4]. Several phase unwrapping algorithms based on branch cuts, selective processing and pixel-wise time-domain algorithms were earlier proposed for improving noise immunity [16–19]. These methods are complicated and have several disadvantages, either they consume large computing time or remain ineffective when large number of phase discontinuities and singularities are present. Noise filtering methods are useful under certain conditions [20] but may remove vital phase gradient information as well.

Here, a novel method of wavefront estimation, which involves a ‘virtual’ Hartmann-Shack (HS) wavefront sensor, is presented. The proposed method attempts to overcome the above mentioned problems by avoiding direct phase unwrapping. The primary advantage of the method lies in the fact that, unlike conventional methods that use wrapped phase images, here, the wavefront reconstruction is performed by using the complex-valued phase factor obtained from the wrapped phase based on detecting the local wavefront slopes with the help of simulated microlenses. The captured interferogram images and wavelength information are sufficient for wavefront reconstruction. There are several parameter tradeoffs in the performance of a conventional HS wavefront sensor [21]. The advantage of the virtual HS unwrapping is that the parameters of the HS can be adjusted without any design difficulties so as to optimise the performance of the digital PS PDI wavefront sensor in real-time. In addition, it is shown that the iterative method converges quickly even in the presence of noise. Related, it was shown earlier that it is possible to obtain unwrapped phase from the derivatives of the wrapped phase despite the phase jumps [22].

It is important to compare the proposed virtual HS based phase unwrapping with existing methods. Since conventional methods [14, 15] failed to be consistent, a Fourier phase unwrapping method [23], which is gaining popularity due to its speed and accuracy benefits and utilizes the complex-valued phase factor for wavefront reconstruction, was implemented on the digital PS PDI interferograms for the first time. The use of least-squares minimization and fastFourier transform (FFT) makes the phase unwrapping method more robust and computationally efficient. A simple and straightforward approach to evaluate the performance of unwrapping is to compare the original wrapped phase with the rewrapped reconstructed wavefront. The rewrapped phase maps obtained from the wavefronts reconstructed using the Fourier-based unwrapping matched well with the raw unwrapped phase maps obtained from the interferograms, giving a clear indication of the effectiveness of the least-squares minimization approach. Another analogous method of obtaining unwrapped phase from digital holograms using Fourier transforms was recently proposed for microscopy applications [24].

In this work, phase-shifted interferograms were captured while sensing five different Zernike aberration modes with the digital PS PDI. They were analysed with the proposed virtual HS and Fourier phase unwrapping methods. Both methods proved to be useful in the context of a digital PS PDI. Through adaptive optics simulation studies in the presence of noise, it is shown that the performance of a virtual HS phase unwrapping method is superior to the Fourier method. Attempts were earlier made to combine a PDI with HS wavefront sensor [25] to utilize the high-resolution capability of PDI and the robust nature of the HS. Through simulations, it is shown that the requirements of robustness and high-resolution can be met by combining Fourier and virtual HS unwrapping methods in a digital PS PDI, thereby potentially reducing the complexity involved in building a dual setup with two different wavefront sensors. The proposed method can be readily adopted by any application involving two-dimensional phase unwrapping.

## 2. Method

#### 2.1. Experimental setup

The expanded beam from a 632.8 nm wavelength He-Ne laser was used as incident light. A MEMS deformable mirror from Boston Micromachines™ with 140 actuators, calibrated with a commercial HS wavefront sensor (Thorlabs WFS150-7AR, lenslet pitch = 150 *μ*m, effective focal length = 5.2 mm, pixel size = 4.65 *μ*m) in closed-loop was used to generate optical aberrations of known magnitude in the aperture plane. The beam diameter in the aperture plane was 4 mm. The beam diameter was reduced to 3 mm on the HS wavefront sensor, equivalently, the wavefront sampling was done using 20 *×* 20 microlenses. Here, the first four orders of Zernike polynomial coefficients estimated by the commercial HS were used. The PS PDI phase mask was addressed on a Holoeye SLM (Model: LC 2012; pixel pitch is 36 *μ*m; pixel fill factor is 58%) and it was placed in the focal plane of a 500 mm lens to mimic a PS PDI. The SLM was calibrated to avoid inherent phase variations by grayscale nonlinearities.

The optical setup is identical to that shown in [4] and an illustration is shown in Fig. 1. The phase mask consists of two squared pinholes at the focal plane: one near the optical axis and another laterally shifted from the center. The fringe visibility and the overall performance of a PDI is affected if the pinhole diameter is larger than the Airy disk diameter when no aberrations are present [26, 27]. Thus, the diameter of the central pinhole was kept lower than 0.19 mm (i.e. lower than 5 SLM pixels). The laterally shifted pinhole, which generates the reference wavefront, needs to be sufficiently small to generate a good quality reference wavefront using a square shaped pinhole. Also, the need for high quality fringes demands an increase in the reference pinhole size to match the diffracting pinhole size. An alternative is to adjust the optical transmission through the SLM by controlling the base grayscale in the two pinholes. However, the use of SLM for phase modulation to achieve PS limits the control on SLM transmission to a certain extent and hence fringe visibility. The fringe visibility is not a fundamental limitation and an SLM with a smaller pixel size can further improve the quality of the reference wavefront and thus the fringes.

Wavefronts were represented using a linear combination of ANSI standard Zernike polynomials (ANSI Z80.28) as follows:

where*a*represents Zernike coefficients and index

_{j}*j*represents the order of Zernike polynomial

*Z*.

_{j}#### 2.2. Wavefront reconstruction with a virtual Hartmann-Shack wavefront sensor

Interferograms obtained with phase shifts of 0, *π*/2, and *π* resulted in three intensity measurements, *I*_{1}, *I*_{2}, and *I*_{3} that may be combined to obtain the wrapped phase as follows [4]:

*ϕ*represents the wrapped phase corresponding to the unwrapped phase

^{w}*ϕ*.

As an example, consider the wrapped phase shown in Fig. 2(b), obtained from Eq. (2), of a wavefront (to be reconstructed) incident on a virtual HS sensor having S*×*S microlenses. The focal length and aperture diameter of the simulated microlenses were chosen to be *f*= 40 mm and *d* = 0.5 mm respectively. The intensity,
${\mathfrak{J}}_{n}$ at the focal plane of the virtual HS corresponding to the *n ^{th}*subaperture can be computed as follows,

*p*(

*ξ, η*)is the subaperture pupil function. ${\varphi}_{n}^{w}$ is the wrapped phase corresponding to the

*n*subaperture. Equation (3) is valid within scalar diffraction theory for monochromatic light. The location of the focal spots was determined using an intensity weighted centroiding algorithm. The local ‘

^{th}*x*’ and ‘

*y*’slopes of the wavefront are obtained using ${s}_{n}^{\overrightarrow{x}}=\mathrm{\Delta}{\overrightarrow{x}}_{n}/f$, where $\mathrm{\Delta}{\overrightarrow{x}}_{n}$ represents the displacement of the focal spots from the centre along ‘x’ and ‘y’ directions. The wavefront is estimated by using the wavefront slope sampling geometry of Southwell [28].

The accuracy of wavefront sensing can be evaluated with the Strehl ratio and root-mean-square (RMS) wavefront error estimated after each wavefront correction. The number of sampling microlenses can be modified for computational efficiency and accuracy by monitoring Strehl ratio and RMS wavefront error. The measurements were repeated on five different daysand the calibration of the deformable mirror with the commercial HS wavefront sensor led to minor changes in the coefficients corresponding to unwanted Zernike aberrations (with a standard deviation of ∼8 nm).

#### 2.3. Fourier phase unwrapping

Fourier-based unwrapping methods are faster and are known to be robust [23, 29]. The Fourier algorithm implemented here finds an optimum unwrapped wavefront in agreement with all observations subject to minimization of the following objective function:

Here, the function *f*(*x,y*) denotes the (unknown) unwrapped wavefront to be recovered,
$\mathrm{\Delta}{\varphi}_{x}\left(x,y\right)={\varphi}^{w}\left(x+p,y\right)-{\varphi}^{w}\left(x,y\right),\mathrm{\Delta}{\varphi}_{y}\left(x,y\right)={\varphi}^{w}\left(x,y+p\right)-{\varphi}^{w}\left(x,y\right)$ and similarly
$\mathrm{\Delta}{f}_{x}\left(x,y\right)=f\left(x+p,y\right)-f\left(x,y\right)$ and
$\mathrm{\Delta}{f}_{y}\left(x,y\right)=f\left(x,y+p\right)-f\left(x,y\right)$. Here, *p* is the sampling period. We will see below that minimization of Eq. (4) is an ill-posed problem in the sense of Hadamard because there is not a unique solution. Even worse, it can yield catastrophic results in the presence of measurement noise. We therefore added a Tikhonov regularization that opts for the solution *f*(*x,y*) with minimum norm. The scalar regularization parameter *γ*has to be chosen with respect to the signal-to-noise ratio (SNR) of the measurements. In all our investigations we used values of γ= 10^{‒12}.

For minimization of Eq. (4) it is convenient to slightly reformulate it by using the convolution operator ⊗ and the odd impulse pairs consisting of Dirac delta functions: *h _{x}*(

*x,y*) =

*δ*(

*x*+

*p,y*) ‒

*δ*(

*x,y*) and

*h*(

_{y}*x,y*) =

*δ*(

*x,y + p*) ‒

*δ*(

*x,y*). The objective function can hence be written:

We then can make use of Parseval’s theorem to look at the problem in Fourier space. When we denote the Fourier transform of a function *g*with *ĝ*we find

*η*are variables in frequency space and

*ĥ*(ν,

_{x}*η*) = i2 sin(πν)exp(iπν) and

*ĥ*(ν,

_{y}*η*) = i2sin(π

*η*)exp(iπ

*η*) are the Fourier transforms of the odd impulse pairs. To find the optimum $\widehat{f}$ that minimizes

*L*we calculate the gradient of

*L*with respect to the function values $\widehat{f}$(ν,

*η*)and set it to zero [8] and,

*ν*and

*η*for the sake of brevity, and ${\widehat{f}}_{R}$ and ${\widehat{f}}_{I}$ are the real and imaginary parts of $\widehat{f}$, respectively. Solving Eq. (7) for $\widehat{f}$(

*ν, η*) is straightforward and yields the Fourier transform of the unwrapped wavefront:

Equation (8) has a pole at *ν* = *η* = 0 that can be avoided by setting δ to values different from zero. The unwrapped wavefront *f*(*x,y*) can be derived from Eq. (8) using the inverse Fourier transform [23].

Finally, it should be noted that in order to calculate reasonable difference values Δ*ϕ _{x}*(

*x,y*) and Δ

*ϕ*(

_{y}*x,y*)along the borders of the rectangular pupil of

*ϕ*(

^{w}*x,y*)we have used the approach ofElster et al. [30], simplifying the problem to a special periodic case (by assuming the wavefront

*f*(

*x,y*) to be inherently periodic).

The resulting wavefront *f*(*x,y*) is potentially different from the initial phase distribution *ϕ ^{w}*(

*x,y*) when wrapped again. The reason is that the differentiation used to obtain Δ

*ϕ*(

_{x}*x,y*) and Δ

*ϕ*(

_{y}*x,y*) can produce inconsistent data if singularities are present on

*ϕ*(

^{w}*x,y*) (for example due to noise, etc.). The differences are usually small. However, in some applications it may not be appropriate for an unwrapping algorithm to change the actual phase data. To avoid this effect, the BIAS-function (27π-steps) of the unwrapped wavefront can be extracted and added to the initial phase values:

Doing so will ensure that the resulting wavefront *ϕ*(*x,y*) equals the measured phase values *ϕ ^{w}*(

*x,y*) when wrapped again.

#### 2.4. Evaluation metrics

Strehl ratio and RMS wavefront error were used to evaluate the performance of phase unwrapping methods. The residual wavefront error (*φ*) is calculated as the difference of the commercial HS measured wavefront and the wavefront reconstructed by the virtual HS and Fourier phase unwrapping methods. Strehl ratio is calculated by taking a ratio of the on-axis intensity at the focal plane due to the complex phase field caused by *φ*and the maximum attainable on-axis intensity at the focal plane due to a flat wavefront. Deviation of the estimated wavefront from the HS measured wavefront would lead to a Strehl ratio lower than 1. The RMS wavefront error is calculated as follows:

It should be noted that in the evaluation of *φ*, tip, tilt and piston terms are not included since they do not degrade the point spread function (PSF) quality. These metrics quantify the accuracy of the digital PS PDI and consequently the quality of the phase unwrapping procedure.

## 3. Results

#### 3.1. Experimental results

Primary defocus, primary astigmatism, coma, secondary astigmatism and a linear combination of defocus and astigmatism were introduced with the deformable mirror. The peak-to-valley of the generated aberrations as measured by the commercial HS wavefront sensor are equal to 1.82 *μ*m, 2.32 *μ*m, 1.85 *μm*, 1.83 *μ*m and 1.86 *μ*minthe same order. Three interferograms (*I*_{1}, *I*_{2} and *I*_{3}) were captured by the CCD camera for each introduced aberration with 0, π/2 and *π* phase shifts introduced by the SLM. The wrapped phase was then evaluated using Eq. (2). The wavefronts were reconstructed using the proposed virtual HS and Fourier unwrapping methods by following the steps illustrated in Sections 2.2 and 2.3 respectively.

The performance of wavefront reconstruction with the proposed virtual HS phase unwrapping and the Fourier phase unwrapping methods in comparison with the wavefronts measured by a commercial HS wavefront sensor is shown in Fig. 2. The wavefronts measured by the commercial HS wavefront sensor are shown in Fig. 2(a). The wrapped phase images, *ϕ ^{w}*(

*x,y*), calculated using Eq. (2), obtained from the digital PS PDI are shown in Fig. 2(b). The wavefronts reconstructed by the virtual HS phase unwrapping are shown in Fig. 2(c). Most ophthalmic and optical microscopy applications deal with Zernike aberrations up to (or below) the fifth order. Hence, the reconstructed wavefronts were Zernike decomposed using least-square fitting [31] with the first 21 Zernike polynomials alone and are shown in Fig. 2(d). Here, tip, tilt and piston terms are excluded. The wavefronts reconstructed using the virtual HS phase unwrapping were rewrapped and are shown in Fig. 2(e). Also, wavefronts reconstructed with Fourier unwrapping method are shown in Fig. 2(f). The corresponding Zernike decomposed wavefronts (first 21 terms excluding tip, tilt and piston) and rewrapped phase images are shown in Figs. 2(g) and 2(h) respectively. It was found that the implemented conventional methods [14,15] were inconsistent while virtual HS and Fourier unwrapping methods were highly repeatable and similar results were obtained for different peak-to-valley of induced aberrations.

Not surprisingly, it was found that the accuracy of wavefront sensing depends on the sampling as shown in Fig. 3. An increase in the accuracy of wavefront sensing is expected with increased wavefront sampling in a HS wavefront sensor [32]. A similar trend was observed in Strehl ratio with increasing wavefront sampling in virtual HS unwrapping as shown in Fig. 3 while sensing defocus, secondary astigmatism and a combination of defocus and astigmatism. However, no significant trend was observed in the case of astigmatism and coma. It can be noted from Fig. 3 that the optimum wavefront sampling, S_{opt}, that results in maximum Strehl ratio is different for the five induced aberrations. Hence the Strehl ratio was estimated for different wavefront sampling (integer steps in the range: 10 ≤S ≤30) for each induced aberration case. The determined values of S_{opt} while sensing defocus, astigmatism, coma, secondary astigmatism and a combination of defocus and astigmatism were found to be 17, 18, 28, 23 and 22 respectively. Also, S_{opt} depends on noise and therefore, it needs to be determined for each iteration. For all practical purposes, where prior information regarding the aberrations is not available, the optical design should incorporate the measurement of Strehl ratio or an equivalent metric to determine optimum wavefront sampling. Since beam resizing does not affect wavefront reconstruction, the aperture diameter (*d* = 0.5 mm) can be fixed while searching for S_{opt}.

The peak-to-valley of the wavefronts estimated by the virtual HS are 2.64 *μ*m, 1.78 *μ*m, 3.05 *μ*m, 1.56 *μ*m and 2.06 *μ*m in the same order. The peak-to-valley of the wavefronts estimated by the Fourier method are 1.86 *μ*m, 1.19 *μ*m, 2.98 *μ*m, 2.91 *μ*m and 1.72 *μ*m. Figure 4 shows the dominant Zernike coefficient values for each Zernike aberration case as estimated by the proposed virtual HS unwrapping in comparison with the Fourier unwrapping method and commercial HS wavefront sensor measurements. Evidently, the virtual HS method is significantly better than Fourier unwrapping when detecting astigmatism and coma. The absolute difference in the estimated defocus coefficient and the commercial HS measured defocus coefficient is only 7% lower in the case of virtual HS as compared to the Fourier method. However, the virtual HS underestimates the secondary astigmatism coefficient to a larger degree than the Fourier method. Consequently, a combination of these two methods while monitoring PSF quality would be helpful in improving convergence in a digital PS PDI.

On average, the estimated Strehl ratio for the five different aberrations sensed using the virtual HS phase unwrapping is 0.09 ± 0.07 and with Fourier unwrapping it was found to be 0.06 ± 0.08. Similarly, the average RMS wavefront error in the case of virtual HS is 0.22 ± 0.07 *μ*m and with Fourier unwrapping it was 0.25 ± 0.08 *μ*m. The virtual HS unwrapping is noticeably better than the Fourier method in terms of both the metrics.

The significant magnitude of the estimated residual aberrations indicate the need for a closed-loop operation as suggested earlier [4]. Through simulation studies, we examine the performance of the virtual HS and Fourier unwrapping methods in the following section.

#### 3.2. Simulation results

Simulations were performed to test the performance of the digital PS PDI. Five different wavefronts (*φ*) including defocus, astigmatism, coma, secondary astigmatism and a combination of defocus and astigmatism were simulated based on Zernike polynomials. The peak-to-valley of the simulated aberrations is the same as that detected by the commercial HS wavefront sensor. The intensity at the detector plane in a digital PS PDI can be expressed as follows:

*P*(

*x, y*) is the circular pupil function and the SLM pinhole mask is defined as follows:

The constant ‘*c*’ helps in adjusting fringe contrast and *θ*_{1}and *θ*_{2}can be controlled to achieve PS. Here, *L* = 108 *μ*m, *l* = 72 *μ*m and *h* = 252 *μ*m, in agreement with the size of the SLM pinholes and spacing between them. *I*_{1}, *I*_{2} and *I*_{3} are evaluated for three phase shifts of 0, π/2 and πusing Eq. (11). Next, the wrapped phase, *ϕ ^{w}*is calculated from Eq. (2). Finally, wavefronts are reconstructed using the virtual HS and Fourier unwrapping methods following the steps described in sections 2.2 and 2.3. Here, S = 20 was chosen for simulation studies since the median of the experimentally derived S

_{opt}values is 20, when coma is excluded. Related, wavefront sampling with the commercial HS is also S = 20.

A simulation-based comparison of the two phase unwrapping methods is shown in Fig. 5. Simulated wavefronts (ϕ) and the wrapped phase images (*ϕ ^{w}*) are shown in Figs. 5(a) and 5(b). The reconstructed wavefronts (after first iteration), Zernike decomposed wavefronts (using Zernike terms 3 ≤

*j*≤ 20) and rewrapped phase images are shown in Figs. 5(c)–Figs. 5(e) for the case of virtual HS phase unwrapping. Similarly, reconstructed wavefronts, Zernike decomposed wavefronts (using Zernike terms 3 ≤

*j*≤ 20) and rewrapped phase images from Fourier unwrapping are shown in Figs. 5(f)–Figs. 5(h). In agreement with experimental results, the average estimated Strehl ratio in the case of virtual HS unwrapping (0.68 ± 0.20) is clearly higher than for the Fourier method (0.63 ± 0.28). Moreover, the RMS wavefront error at the end of first iteration for virtual HS unwrapping is 0.18 ± 0.04

*μ*m and for the Fourier method it is 0.17 ± 0.03

*μ*m. It can be noted that the Strehl ratio values obtained through simulations are higher in comparison with the experimental values. This is attributed to the absence of noise in the simulated interferograms.

In most practical situations, optical aberrations are expressed as a linear combination of Zernike polynomials and to study such cases, iterative calculations were performed on 100 randomly simulated wavefronts and the residual wavefront error, *φ* obtained at the end of a given iteration was assumed to be the new wavefront, ϕfor the following iteration. Random wavefronts were generated using Eq. (1) and the Zernike coefficients, *a _{j}* (

*j*≤ 14) are randomly picked from a uniform distribution with mean 0 and variance 1. Here too, piston, tip and tilt coefficients were assumed to be equal to zero. Strehl ratio and RMS wavefront error were estimated after each iteration. Figure 6 shows the effect of wavefront sampling in a virtual HS unwrapping in terms of Strehl ratio and RMS wavefront error. Clearly, an increase in wavefront sampling increased the accuracy of wavefront sensing and the difference between S = 5 and S = 10 cases is prominent in the first few iterations and convergence was observed beyond four iterations. Increasing the wavefront sampling beyond S = 10 did not show any significant increase in Strehl ratio. This is due to the fact that only aberrations up to the fourth order of Zernike are introduced here. When high-order aberrations beyond the fourth order were included, the optimal wavefront sampling was found to be higher. However, wavefront sampling cannot be indiscriminately increased for two reasons. Firstly, due to increased computing time and secondly, in the presence of noise, higher-order Zernike coefficients erroneously assume large values.

A comparison of the Strehl ratio values for virtual HS and Fourier unwrapping methods is shown in Fig. 7(a). The calculated Strehl ratio at the end of first iteration in the case of Fourier unwrapping is higher than the virtual HS in the case of 100 randomly simulated aberrations. This can be attributed to the presence of high spatial frequency components in the aberrated phase maps. However, the Strehl ratio converges to a lower value in the Fourier case in comparison with the virtual HS unwrapping. Given that the convergence takes a different path in both these methods, it would be interesting to see the convergence path when both the methods are combined in such a way that the method which yields a better Strehl ratio is chosen for each step of the correction. Figure 7(a) suggests that combining the methods helps in increasing the convergence efficiency.

To study the effect of noise, white Gaussian noise was added to the simulated interferograms (*I*_{1}, *I*_{2} and *I*_{3}). A comparison of the Strehl ratio for three different cases of SNR 10, 20 and 30 dB is shown in Fig. 7(b). The performance of virtual HS is clearly superior to Fourier unwrapping in the presence of noise. When the SNR is 20 dB and lower, the Fourier method fails to reach the Marechal’s criterion whereas diffraction-limited performance is achieved in the virtual HS phase unwrapping method at the end of five iterations although the SNR is as low as 10 dB. From this result, it may be noted that under the situations of high noise, combining these two methods will not improve convergence since the calculations made with the help of Fourier unwrapping method will remain redundant.

## 4. Discussion

In comparison with the wavefronts sensed by a commercial HS wavefront sensor shown in Fig. 2(a), the reconstructed wavefronts in both the implemented phase unwrapping techniques shown in Figs. 2(c) and Figs. 2(f) seem to be questionable. It has to be noted that the measurement noise allows the appearance of phase singularities [33] in the wrapped phase maps as clearly visible in Fig. 2(b), which make the result of phase unwrapping ambiguous [34]. In addition, inaccuracies could be caused by the generation of the reference wavefront with a finite-sized square pinhole on the SLM, unequal tilt aberrations due to non-common path measurements by the commercial HS and digital PS PDI wavefront sensors, and the presence of noise in the interferograms. Furthermore, in the presence of noise, wavefront slopes are known to be underestimated [20]. However, after fitting the reconstructed wavefronts using the first 21 Zernike polynomials, the resultant decomposed wavefronts match reasonably well the aberrations sensed by the commercial HS. The use of median filtering on the experimentally obtained interferograms did not show a statistically significant increase in the accuracy of wavefront sensing. Alternative filtering methods from the literature which were not implemented here may improve the accuracy. Yet, utilizing all the information that the interferograms provide will reduce errors due to unwanted smoothing. Moreover, choosing a lower wavefront sampling in a virtual HS method is an indirect mechanism of filtering the high spatial frequencies in the reconstructed phase maps. This is evident when comparing the reconstructed wavefronts in Figs. 2(c) and Figs. 2(f). The Fourier unwrapping method left behind traces of the locations of phase discontinuities in the wrapped phase. However, this was not prominent in the virtual HS case. This is because in the Fourier case, the reconstruction is based on all observations in the wrapped phase and in the virtual HS case, the sampling is generally much lower than the number of active pixels in the interferograms and hence these features are filtered out. In addition, the proposed method includes Zernike decomposition, which helps in representing the reconstructed wavefront in terms of lower-order Zernike polynomials alone, functioning like a Zernike filter.

Here, we made use of the intensity weighted centroiding method, which is a computationally efficient and accurate focal spot centroiding algorithm. However, it should be noted that since this centroiding method detects the location of peak intensity rather than the center-of-mass, the wavefront sampling should be cleverly chosen such that the phase corresponding to individual subapertures contain predominantly tilt aberration alone. If the local wavefront has aberrations other than tilt, the simulated focal spots will not be diffraction-limited and therefore, the intensity weighted centroiding algorithm may fail. Although the use of advanced centroiding methods might increase computational time, they could be reliable while changing slope sampling. However, the use of matched filter and iteratively weighted center-of-gravity methods [21, 35] in virtual HS did not improve the accuracy of wavefront sensing when applied on experimentally obtained interferograms. The accuracy of centroiding improves with increasing iteration number, which is attributed to the reduced magnitude of aberrations.

The lower limit to wavefront sampling in virtual HS unwrapping is set by the Nyquist sampling criterion and depends on the number of Zernike modes of interest. The upper limit is set by the number of pixels in the interferograms, computational constraints and noise. While determining the optimum wavefront sampling, only the relative quality of wavefront correction is of interest. Hence, wavefront sampling can be optimized in a practical situation by taking light power measurements behind a pinhole [2], whose diameter is nearly equal to the Airy disk diameter [36]. Since the accuracy of wavefront sensing is dependent on wavefront sampling, rewrapping of the reconstructed wavefront depends on wavefront sampling too. Unlike Fourier unwrapping method where every pixel in the interferograms is used for residual error minimization, the inherent spatial filtering due to wavefront sampling in a virtual HS unwrapping does not allow rewrapping to be a good measure of the wavefront sensing accuracy.

Due to alignment and imperfections in closed-loop calibration of the deformable mirror with the commercial HS wavefront sensor, it is not possible to introduce pure Zernike aberrations with the deformable mirror and a few lower-order Zernike coefficients assume significant nonzero values. For instance, while generating primary defocus, the RMS wavefront error due to the unwanted Zernike aberrations was measured to be ~10% of the total RMS wavefront error. Deviations in the estimated Zernike coefficients from the commercial HS measured values could be caused to a certain degree by differences in contribution from other Zernike aberrations. The use of an SLM with greater fill-factor and smaller pixel size could improve tunability and accuracy of the generated reference wavefront [4]. The number of phase-shifting steps could be increased to partly reduce the effects of noise on wavefront sensing at the cost of increased data acquisition time.

The increase in Strehl ratio with increasing wavefront sampling shown in Fig. 3 was not observed in the case of primary astigmatism and coma aberrations. Also, for the aberrations shown in Fig. 3, for S > 25, a drop in the Strehl ratio was observed. This is attributed to the presence of noise. The resulting phase unwrapping errors lead to unexpected oscillations in Fig. 3. However, as expected, the plots of RMS wavefront error corresponding to the Strehl ratio plots in Fig. 3 showed a decreasing trend. Since aberrations and noise affect the determination of optimal sampling, it is recommended to evaluate S_{opt} for every single iteration. The reliability of the estimated optimal sampling improves with increasing Strehl ratio and iteration number.

The use of a lower number of pixels per subaperture in a virtual HS unwrapping method can improve speed. However, in simulation studies, a decrease in the number of pixels per subaperture from 32 to 16 showed that the mean Strehl ratio for 100 random wavefronts after the first iteration drops from 0.64 to 0.59. Our simulations excluded possible optical misalignments and wavefront compensation errors. As can be seen from Fig. 7, in the presence of noise, virtual HS dominates Fourier unwrapping. Experimental inferences suggest that combining the two methods improves wavefront sensing accuracy in certain cases, unless it is verified that one of the methods is entirely redundant. Simulation results showed that the presence of noise makes the virtual HS method more robust than the Fourier unwrapping method.

Although the digital PS PDI was used as an example to test the working of the proposed virtual HS phase unwrapping method, it can potentially be applied to any application involving phase unwrapping (with or without noise) including interferometry. The novel virtual HS unwrapping method is worth applying in situations where the incoming aberrations are dominated by low-order Zernike polynomials in the likes of vision science and microscopy. For applications requiring high precision phase unwrapping, a combination of virtual HS and Fourier unwrapping will be helpful in achieving good accuracies.

The unique advantage of the virtuality is that it allows to alter several HS wavefront sensor parameters so as to optimize wavefront reconstruction for a given aberration and noise. In addition, it provides an alternative to phase unwrapping and the wavefront sensing results of virtual HS and Fourier phase unwrapping methods may be combined to further optimize performance of the digital PS PDI. The most important parameter is the wavefront sampling, which can be increased to sense the wavefront more closely although ultimately this will be limited by measurement noise and computational time for wavefront reconstruction. It has to be noted that the proposed technique could cause errors due to discretization arising from the sampling of the wavefront slopes. The choice of wavefront sampling can be varied based on the application at hand.

In conclusion, a novel method of phase unwrapping has been proposed using the concept of a virtual HS, which uses the complex phase function to determine the local wavefront slopes. The method has been compared with a robust Fourier unwrapping method. It has been shown through simulations that the proposed method is superior in comparison with the Fourier method in the presence of noise.

## Acknowledgments

The research leading to these results has received funding from Consejería de Educación, Juventud y Deporte of Comunidad de Madrid and the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme ( FP7/2007-2013) under REA grant agreement n°291820. This study was supported by Spanish Government grant FIS2011-25637 to Prof. Susana Marcos. This work was also supported by Science Foundation Ireland (grants: 07/SK/B1239a and 08/IN.1/B2053) and University College Dublin (seed funding: SF665).

## References and links

**1. **V. Akondi, S. Castillo, and B. Vohnsen, “Digital pyramid wavefront sensor with tunable modulation,” Opt. Express **21**(15), 18261–18272 (2013). [CrossRef] [PubMed]

**2. **A. R. Jewel, V. Akondi, and B. Vohnsen, “A direct comparison between a MEMS deformable mirror and a liquid crystal spatial light modulator in signal-based wavefront sensing,” J. Eur. Opt. Soc. Rapid Pub. **8**, 13073 (2013). [CrossRef]

**3. **V. Akondi, S. Castillo, and B. Vohnsen, “Multi-faceted digital pyramid wavefront sensor,” Opt. Commun. **323**, 77–86 (2014). [CrossRef]

**4. **V. Akondi, A. R. Jewel, and B. Vohnsen, “Digital phase-shifting point diffraction interferometer,” Opt. Lett. **39**(6), 1641–1644 (2014). [CrossRef] [PubMed]

**5. **V. Akondi, A. R. Jewel, and B. Vohnsen, “Closed-loop adaptive optics using a spatial light modulator for sensing and compensating of optical aberrations in ophthalmic applications,” J. Biomed. Opt. **19**(9), 096014 (2014). [CrossRef]

**6. **C. Falldorf, M. Agour, C. v. Kopylow, and R. B. Bergmann, “Phase retrieval by means of a spatial light modulator in the Fourier domain of an imaging system,” Appl. Opt. **49**(10), 1826–1830 (2010). [CrossRef] [PubMed]

**7. **C. Falldorf, C. v. Kopylow, and R. B. Bergmann, “Wave field sensing by means of computational shear interferometry,” J. Opt. Soc. Am. A **30**(10), 1905–1912 (2013). [CrossRef]

**8. **U. Schnars, C. Falldorf, J. Watson, and W. Jüptner, *Digital Holography and Wavefront Sensing* (Springer, 2015).

**9. **P. Hariharan, “Digital phase-stepping interferometry: effects of multiply reflected beams,” Appl. Opt. **26**(13), 2506–2507 (1987). [CrossRef] [PubMed]

**10. **O. Y. Kwon, “Multichannel phase-shifted interferometer,” Opt. Lett. **9**(2), 59–61 (1984). [CrossRef] [PubMed]

**11. **H. Kadono, N. Takai, and T. Asakura, “New common-path phase shifting interferometer using a polarization technique,” Appl. Opt. **26**(5), 898–904 (1987). [CrossRef] [PubMed]

**12. **C. R. Mercer and K. Creath, “Liquid-crystal point-diffraction interferometer for wave-front measurements,” Appl. Opt. **35**(10), 1633–1642 (1996). [CrossRef] [PubMed]

**13. **S. Heshmat, S. Tomioka, and S. Nishiyama, “Performance evaluation of phase unwrapping algorithms for noisy phase measurements,” Int. J. Optomechatronics , **8**(4), 260–274 (2014). [CrossRef]

**14. **M. A. Herráez, D. R. Burton, M. J. Lalor, and D. B. Clegg, “Robust, simple, and fast algorithm for phase unwrapping,” Appl. Opt. **35**(29), 5847–5852 (1996). [CrossRef] [PubMed]

**15. **R. M. Goldstein, H. A. Zebker, and C. L. Werner, “Satellite radar interferometry: two-dimensional phase unwrapping,” Radio Sci. **23**(4), 713–720 (1988). [CrossRef]

**16. **J. M. Huntley, “Noise-immune phase unwrapping algorithm,” Appl. Opt. **28**(16), 3268–3270 (1989). [CrossRef] [PubMed]

**17. **J. M. Huntley and H. Saldner, “Temporal phase unwrapping algorithm for automated interferogram analysis,” Appl. Opt. **32**(17), 3047–3052 (1993). [CrossRef] [PubMed]

**18. **J. A. Quiroga and E. Bernabeu, “Phase-unwrapping algorithm for noisy phase-map processing,” Appl. Opt. **33**(29), 6725–6731 (1994). [CrossRef] [PubMed]

**19. **R. Cusack, J. M. Huntley, and H. T. Goldrein, “Improved noise-immune phase-unwrapping algorithm,” Appl. Opt. **34**(5), 781–789 (1995). [CrossRef] [PubMed]

**20. **J-S. Lee, K. P. Papathanassiou, T. L. Ainsworth, M. R. Grunes, and A. Reigber, “A new technique for noise filtering of SAR interferometric phase images,” IEEE Trans. Geosci. Remote Sens. **36**(5), 1456–1465 (1998). [CrossRef]

**21. **V. Akondi and B. Vohnsen, “Myopic aberrations: impact of centroiding noise in Hartmann Shack wavefront sensing,” Ophthalmic Physiol. Opt. **33**(4), 434–443 (2013). [CrossRef] [PubMed]

**22. **J Arines, “Least-squares modal estimation of wrapped phases: application to phase unwrapping,” Appl. Opt. **42**(17), 3373–3378 (2003). [CrossRef] [PubMed]

**23. **M. D. Pritt and J. S. Shipman, “Least-squares two-dimensional phase unwrapping using FFT’s,” IEEE Trans. Geosci. Remote Sens **32**(3), 706–708 (1994). [CrossRef]

**24. **I. Iglesias, “Phase estimation from digital holograms without unwrapping,” Opt. Express **22**(18), 21340–21346 (2014). [CrossRef] [PubMed]

**25. **J. M. Bueno, E. Acosta, C. Schwarz, and P. Artal, “Wavefront measurements of phase plates combining a point-diffraction interferometer and a Hartmann-Shack sensor,” Appl. Opt. **49**(3), 450–456 (2010). [CrossRef] [PubMed]

**26. **R. N. Smartt and W. H. Steel, “Theory and application of point-diffraction interferometers,” J. Appl. Phys. **14**(141), 351–356 (1975). [CrossRef]

**27. **F. Bai, X. Wang, K. Huang, N. Gu, S. Gan, and F. Tian, “Analysis of spatial resolution and pinhole size for single-shot point-diffraction interferometer using in closed-loop adaptive optics,” Opt. Commun. **297**, 27–31 (2013). [CrossRef]

**28. **W. H. Southwell, “Wave-front estimation from wave-front slope measurements,” J. Opt. Soc. Am. **70**(8), 998–1006 (1980). [CrossRef]

**29. **V. V. Volkov and Y. Zhu, “Deterministic phase unwrapping in the presence of noise,” Opt. Lett. **28**(22), 2156–2158 (2003). [CrossRef] [PubMed]

**30. **C. Elster and I. Weingärtner, “Solution to the shearing problem,” Appl. Opt. **38**(23), 5024–5031 (1999). [CrossRef]

**31. **C. C. Paige and M. A. Saunders, “LSQR: An algorithm for sparse linear equations And sparse least squares,” ACM Trans. Math. Soft. **8**, 43–71 (1982). [CrossRef]

**32. **R. M. Basavaraju, V. Akondi, S. J. Weddell, and R. P. Budihal, “Myopic aberrations: Simulation based comparison of curvature and Hartmann Shack wavefront sensors,” Opt. Commun. **312**, 23–30 (2014). [CrossRef]

**33. **N. B. Baranova, A. V. Mamaev, N. F. Pilipetsky, V. V. Shkunov, and B. Y. Zel’dovich, “Wave-front dislocations: topological limitations for adaptive systems with phase conjugation,” J. Opt. Soc. Am. **73**(5), 525–528 (1983). [CrossRef]

**34. **J. M. Huntley, “Three-dimensional noise-immune phase unwrapping algorithm,” Appl. Opt. **40**(23), 3901–3908 (2001). [CrossRef]

**35. **V. Akondi, R. M. Basavaraju, and R. P. Budihal, “Centroid detection by Gaussian pattern matching in adaptive optics,” Int. J. Comp. **26**(1), 30–35 (2010).

**36. **A. R. Jewel, V. Akondi, and B. Vohnsen, “Optimization of sensing parameters for a confocal signal-based wave-front corrector in microscopy,” J. Mod. Opt. **62**(10), 786–792 (2015). [CrossRef]