Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group
Journal Home About Issues in Progress Current Issue All Issues Early Posting

Phase unwrapping algorithm based on phase diversity wavefront reconstruction and virtual Hartmann–Shack technology

Ying Zhang, Hua Bao, Naiting Gu, Shuqi Li, Yiqun Zhang, and Changhui Rao

Open Access Open Access

Abstract

Phase unwrapping (PU) algorithms play a crucial role in various phase measurement techniques. Traditional algorithms cannot work well in strong noise environments, which makes it very difficult to obtain the accurate absolute phase from the noisy wrapped phase. In this Letter, we introduce a novel, to the best of our knowledge, phase unwrapping algorithm named PD-VHS. This algorithm innovatively employs point spread function (PSF) filtering to eliminate noise from the wrapped phase. Furthermore, it combines a phase diversity (PD) wavefront reconstruction technology with a virtual Hartmann–Shack (VHS) technology for phase reconstruction and phase unwrapping of the filtered PSFs. In simulations, hundreds of random noise wrapped phases, containing the first 45 Zernike polynomials (excluding piston and the two tilt terms) and the wavefront RMS = 0.5λ and 1λ, are used to compare the classical quality-map guided algorithm, the VHS algorithm with decent noise immunity, with our PD-VHS algorithm. When signal-to-noise ratio (SNR) drops to just 2 dB, the mean root mean square errors (RMSEs) of the residual wavefront between the unwrapped result and the absolute phase of the quality-map guided algorithm and the VHS algorithm are up to 3.99λ, 0.44λ, 4.29λ, and 0.85λ, respectively; however, our algorithm RMSEs are low: 0.11λ and 0.17λ. Simulation results demonstrated that the PD-VHS algorithm significantly outperforms the quality-map guided algorithm and the VHS algorithm under large-scale noise conditions.

© 2024 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

Introduction. The wrapped phase confines absolute phase values within [-$\pi$, $\pi$] and is prevalent in phase measurement-based techniques such as optical metrology [1], digital holography [2], and synthetic aperture radar (SAR) interferometry [3]. Phase unwrapping (PU) is critical for retrieving the absolute phase. The relationship between the absolute and wrapped phases is defined by Eq. (1):

$$\varphi(x,y) = \phi(x,y) + 2\pi k(x,y) + N(x,y),$$
where $\phi$ is the absolute phase, $\varphi$ is the wrapped phase, $(x,y)$ denotes the pixel coordinates, $k$ is the wrap count, and $N(x,y)$ is the noise. The absolute phase $\phi (x,y)$ often models as a sum of Zernike polynomials:
$$\phi(x,y) = \sum_{i=1}a_{i}Z_i(x,y),$$
where $i$ is the order of Zernike polynomial $Z_i$ and $a_{i}$ is Zernike coefficients.

In theory, absent noise or discontinuities, one can accurately recover an absolute phase from $\varphi$. However, noise complicates this process, necessitating robust PU algorithms. These include traditional spatial methods, such as path-tracing using the Itoh condition, with examples such as Goldstein’s branch cut [4], quality-map guided algorithms [5–7], and the mask-cutting algorithm [8]. These algorithms focus on selecting optimal integration paths for accurate PU, a task that becomes difficult in strong noisy conditions. Alternatively, path-independent methods [9,10] minimize phase gradient discrepancies through global optimization, effectively ensuring phase continuity but risking error propagation. While these traditional methods are effective against small-scale noise, their performance declines as noise increases.

Recent advancements in deep learning-based PU methods address PU problem under middle-scale noise conditions, leveraging end-to-end learning models [11–13]. However, these methods need extensive training datasets, which is challenging when data is scarce or costly, and may not perform well under unrepresented conditions.

Distinctly, the “virtual” Hartmann–Shack (VHS) algorithm presented in Ref. [14] enhances noise immunity via wavefront sampling and Zernike mode reconstruction. As depicted in Fig. 1, the VHS algorithm effectively unwraps phases under small-scale to middle-scale noise conditions [Figs. 1(b) and 1(e)]. However, under large-scale noise conditions, spot diffusion [Fig. 1(h)] will result in significant errors in centroid calculation and decrease in PU accuracy. Given these challenges, developing more robust PU methods is essential for applications with large-scale noise, such as magnetic resonance imaging (MRI) [15] where the SNR may be below 5 dB.

 figure: Fig. 1.

Fig. 1. Wrapped phases in the presence of different scales of noise and their corresponding HS spot arrays and PSFs. (a)–(c) Corresponding to the small scale, (d)–(f) corresponding to the middle scale, and (g)–(i) corresponding to the large scale.

Download Full Size | PDF

Here, a novel PU algorithm, denoted as PD-VHS, integrating phase diversity (PD) wavefront reconstruction with the VHS technology, effectively managing PU problems under large-scale noise conditions, is proposed. The innovations of this algorithm lie in the following two points. The primary innovation is its novel approach to filter out noise from the wrapped phase using point spread function(PSF), which provides the possibility for the VHS algorithm to perform PU under large-scale noise conditions. The adoption of this new filtering approach stems from the following reasons: Due to the operation of wavefront sampling, the signal-to-noise ratio (SNR) of the wavefront within a single sub-aperture will be much lower than the SNR of the entire wavefront. Therefore, under large-scale noise conditions, the VHS spot array disperses, while the far-field PSF [Fig. 1(i)] still retains the essential spot morphology information just as the PSFs [Figs. 1(c) and 1(f)] under small-scale and middle-scale noise conditions. However, the presence of a large-scale noise in the wrapped phase introduces high spatial frequency noisy components in the PSF [Fig. 1(i)], hindering effective utilization of spot morphology information. To mitigate this, we have developed a PSF filtering approach, detailed in Sect. $2$, which effectively reduces high-frequency noise while preserving the essential spot details. The secondary innovation of this algorithm is that we use the PD wavefront reconstruction technology [16–18] to reconstruct the wrapped phase with the filtered PSFs and reconstruct the absolute phase with VHS. The reason why we choose PD is that using a single-filtered PSF for phase reconstruction is faced with multi-solution problems and the filtered far-field may not correspond to the original near-field any more, whereas PD can constrain the solution space to ensure the correctness of the solution. Specifically, we introduce defocus aberration to generate modulated PSFs and employ the Gerchberg–Saxton (GS) algorithm for phase reconstruction [19].

To verify and compare the performance of the PD-VHS algorithm, we conducted PU simulations under different noise scales using the traditional quality-map guided, VHS and our PD-VHS algorithms. The implementation of the quality-map guided algorithm follows the principles outlined in Ref. [5]. The VHS algorithm and the PD-VHS algorithm are implemented in accordance with the principles detailed in the Method section. Specifically, the root mean square error (RMSE) of the residual wavefront between the unwrapped phase and the corresponding absolute phase is adopted to evaluate the performance of these algorithms. Simulation results demonstrate that under large-scale noise conditions, the mean RMSE of the PD-VHS algorithm is much lower than the quality map-guided algorithm and the VHS algorithm. Meanwhile, the PD-VHS algorithm is more stable and more accurate than the other two algorithms when the scale of noise is not too large.

The Method section and the Simulations section of this Letter will detail the PD-VHS algorithm’s principles and simulations, respectively, including parameter settings and performance comparisons. The Conclusion section summarizes the findings and implications of this Letter.

Method. The PD-VHS algorithm effectively handles PU under large-scale noise conditions through three key techniques: PSF filtering, PD, and VHS. This algorithm often requires multiple iterations to achieve accurate PU due to imperfections in PSF filtering and errors in PD phase reconstruction. The specific principle of the PD-VHS algorithm is illustrated in Fig. 2.

 figure: Fig. 2.

Fig. 2. Principle of the PD-VHS algorithm.

Download Full Size | PDF

In Fig. 2, Phase O, the object phase to be unwrapped, is the algorithm’s input. In Step 1, a defocus aberration is added to the wrapped phase to create the “modulated wrapped phase.” The corresponding PSFs are then generated according to the scalar diffraction theory, as described in Eq. (3):

$$PSF=|F(Pupil(x,y)\cdot e^{\left(j\cdot \left( \phi(x,y)\right)\right)})|^2,$$
where $Pupil(x,y)$ represents the pupil function and $F$ represents the Fourier transform. In Step 2, we apply filtering techniques to the PSFs, involving threshold reduction [20] and a Butterworth low-pass filter. The PSF to be filtered, as depicted in Fig. 3, is processed according to Eqs. (4), (5):
$$\ PSF=Butterworth\_filter[PSF-T]\quad((PSF < 0) = 0)$$
$$T = \alpha \times mean2(A+B+C+D) + \beta \times std2(A+B+C+D).$$
$A$, $B$, $C$, and $D$ share the same dimensions, and their sizes, along with the parameters $\alpha$ and $\beta$, are adjustable based on the specific requirements of different tasks. For our application, simulations have demonstrated that optimal performance is achieved when the dimensions of $A$ are $1/4$ those of the PSF image, with $\alpha$ set at 2 and $\beta$ at 4. Following threshold reduction, which reduces but does not eliminate noise in the PSF, a Butterworth low-pass filter is applied to further attenuate high-frequency noise while preserving essential image details [21]. In Step 3, the GS algorithm is used for wrapped phase reconstruction with the filtered PSFs, employing both the phase difference convergence criterion and a maximum number of iterations to ensure robust convergence. The implementation details are provided in the pseudocode outlined in Algorithm 1:

Tables Icon

Algorithm 1. GS algorithm

View Table
 figure: Fig. 3.

Fig. 3. Sample of PSF to be filtered.

Download Full Size | PDF

In Step 4, we employ the VHS algorithm based on the Zernike mode to unwrap the wrapped phase obtained from the GS algorithm. Its principle can be expressed as in Eqs. (6), (7):

$$A=D^{+} G$$
$$\left[\begin{array}{c} a_1 \\ a_2 \\ \vdots \\ a_v \end{array}\right]=\left[\begin{array}{cccc} Z_{11} & Z_{12} & \cdots & Z_{1 v} \\ Z_{21} & Z_{22} & \cdots & Z_{2 v} \\ \vdots & \vdots & \ddots & \vdots \\ Z_{2 u, 1} & Z_{2 u, 2} & \cdots & Z_{2 u, v} \end{array}\right]^{+}\left[\begin{array}{c} g_1 \\ g_2 \\ \vdots \\ g_{2 u} \end{array}\right],$$
where $A$ is a column vector composed of Zernike coefficients and $D$ is the recovery matrix, derived from the structural parameters of the HS wavefront sensor. The matrix $D^{+}$ represents the generalized inverse of $D$ and $G$ is a column vector detailing the slopes in the $XY$ direction for each sub-aperture. The variables $u$ and $v$ denote the total numbers of sub-apertures and Zernike coefficients to be recovered, respectively. The slopes within the sub-aperture are calculated through the traditional center of gravity (CoG) method [22]. Step 5 involves accumulating the unwrapped phase from each iteration to form “Phase A,” starting from an initial phase of zero. This process checks whether the specified number of iterations has been reached. If so, the computation halts and outputs the unwrapped phase “Phase A.” If not, the algorithm subtracts “Phase A” from the original wrapped phase “Phase O” and repeats the above steps as necessary.

Simulations. For an optical system, the full width at half maximum (FWHM) of the far-field PSF under parallel light serves as an indicator of the spot sampling rate and is calculated by Eq. (8):

$$\mathrm{FWHM}=\frac{\lambda \cdot \mathrm{f}}{\mathrm{D} \cdot \text{ pixel, }}$$
where $\lambda$, $f$, $D$, and $pixel$ are the wavelength of incident light, the focal length, the diameter of the lens, and the size of the camera pixel, respectively. Here, we use the FWHM to characterize the optical system parameters in the simulations directly. In the simulated optical system for PD wavefront reconstruction, the FWHM is of 2 pixels. The peak to valley (PV) of the defocus aberration is $1\lambda$($1\lambda = 632.8$ nm) and the sizes of the phase map and PSF are $256 \times 256$ and $512 \times 512$, respectively. In the simulated optical system for VHS, the FWHM of each microlens is of 3.5 pixels.

In our simulations, we generate wrapped phases through the following steps: Initially, we generate the first 45 Zernike coefficients, excluding the piston and the two tilt terms, according to the Kolmogorov turbulence theory [23]. We then derive the corresponding absolute phase using Eq. (2) and normalize its RMS to either $0.5 \lambda$ or $1 \lambda$. The phase is then wrapped according to Eq. (1), and various scales of the white Gaussian noise are added to produce noisy wrapped phases. To evaluate and compare the performance of different algorithms under varying scales of aberrations and noise conditions, the quality-map guided algorithm, VHS algorithm, and PD-VHS algorithm were employed to perform PU on these phases. Since the wrapped phase includes only the first 4–45 Zernike polynomials, we specifically restore these coefficients in both the VHS and PD-VHS algorithms. Additionally, given that only the 4th to 45th Zernike terms need to be reconstructed, we selected a $16\times 16$ wavefront sampling rate in both the VHS and PD-VHS algorithms to sample the wrapped wavefront, which is sufficient for accurate phase reconstruction of wrapped phases containing the first 4–45 Zernike polynomials. The PU results are depicted in Fig. 4.

 figure: Fig. 4.

Fig. 4. Simulations: evaluation of quality-map guided, VHS, and PD-VHS algorithms for two scales of (a) $0.5 \lambda$ and (b) $1 \lambda$ aberrations with different scales of noise (SNR is shown in decibels).

Download Full Size | PDF

Obviously, when the SNR is higher than 10 dB, three PU algorithms can unwrap the wrapped phase effectively. However, when the SNR is lower than 10 dB and higher than 2 dB, the quality-map guided algorithm has failed, and the performance of the VHS algorithm decreased dramatically, while the PD-VHS algorithm can still achieve high PU accuracy for different scales of aberrations. When the SNR drops to 2 dB, for the two scales of the absolute wavefront, the mean RMSEs of the PD-VHS algorithm can still achieve the accuracy of $0.11\lambda$ and $0.17\lambda$, while the mean RMSEs of the quality-map guided algorithm and the VHS algorithm are up to $3.99\lambda$ and $0.44\lambda$ and $4.29\lambda$ and $0.85\lambda$. A PU example in this case is shown in Fig. 5.

 figure: Fig. 5.

Fig. 5. Example of PU results of quality-map guided, VHS, and PD-VHS algorithms. (a) Absolute phase (RMS = $1 \lambda$), (b) wrapped phase, (c)–(e) PU results of quality-map guided (RMSE = $4.91\lambda$), VHS (RMSE = $0.89\lambda$), and PD-VHS (RMSE = $0.07\lambda$) algorithms, (f)–(j) PSFs corresponding to phases (a)–(e).

Download Full Size | PDF

These results demonstrated that the PD-VHS algorithm significantly outperforms the quality-map guided algorithm and the VHS algorithm under larger-scale noise conditions. However, the performance of the PD-VHS algorithm has decreased slightly when the SNR drops to 2 dB. That is because the ability of the PD technology to reconstruct residual aberrations is limited and unstable when the noise is too strong. So better PSF filtering methods need to be explored to further improve the accuracy of this algorithm.

Conclusion. In this Letter, we have introduced a novel PU algorithm, termed the PD-VHS algorithm, which innovatively utilizes PSFs to filter noise from wrapped phases. This method addresses PU challenges under large-scale noise conditions by combining PD with VHS wavefront reconstruction technologies. Simulation results reveal that the PD-VHS algorithm significantly surpasses the quality-map guided and VHS algorithms in both PU accuracy and noise immunity. Compared with the deep learning-based PU methods, the PD-VHS algorithm not only avoids the limitations discussed in the introduction, such as dependency on large datasets and potential generalization issues, but also provides a transparent, principle-driven alternative to the “black box” nature of deep learning models, enhancing troubleshooting, interpretation, and implementation ease. Besides, the PD-VHS algorithm does not require a complex network setup or significant computational resources, making it more accessible and easier to integrate into various application scenarios. However, there are areas where the PD-VHS algorithm could be improved. Primarily, its performance degrades under conditions of intense noise, suggesting a need for more effective PSF filtering methods. Additionally, due to its reliance on the GS algorithm and iterative calculations, the process can be time-consuming. Optimizing the algorithm’s structure to reduce computational time is necessary to enhance its practical applicability and efficiency.

Funding

Laboratory Innovation Foundation of the Chinese Academy of Science (YJ22K002); National Natural Science Foundation of China (11727805).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available but may be obtained from the authors upon reasonable request.

REFERENCES

1. T. Zhang, S. Jiang, Z. Zhao, et al., Opt. Express 27, 23173 (2019). [CrossRef]  

2. G. Pedrini, I. Alexeenko, W. Osten, et al., Appl. Opt. 42, 5846 (2003). [CrossRef]  

3. L. Zhou, H. Yu, V. Pascazio, et al., IEEE Trans. Geosci. Remote Sensing 60, 5221510 (2022). [CrossRef]  

4. J. M. Bioucas-Dias and G. Valadao, IEEE Trans. on Image Process. 16, 698 (2007). [CrossRef]  

5. M. A. Herráez, D. R. Burton, M. J. Lalor, et al., Appl. Opt. 41, 7437 (2002). [CrossRef]  

6. M. Arevalillo-Herráez, F. R. Villatoro, and M. A. Gdeisat, IEEE Trans. on Image Process. 25, 2601 (2016). [CrossRef]  

7. X. Long, H. Bao, C. H. Rao, et al., Opto-Electron. Eng. 47, 200111 (2020). [CrossRef]  

8. X. Su and W. Chen, Opt. Lasers Eng. 42, 245 (2004). [CrossRef]  

9. J. Strand, T. Taxt, and A. K. Jain, IEEE Trans. on Image Process. 8, 375 (1999). [CrossRef]  

10. M. D. Pritt and J. S. Shipman, IEEE Trans. Geosci. Remote Sensing 32, 706 (1994). [CrossRef]  

11. G. Dardikman-Yoffe, D. Roitshtain, S. K. Mirsky, et al., Biomed. Opt. Express 11, 1107 (2020). [CrossRef]  

12. J. Liu, Q. Wu, X. Sui, et al., PhotoniX 2, 5 (2021). [CrossRef]  

13. J. Chen, Y. Kong, D. Zhang, et al., Opt. Express 31, 29792 (2023). [CrossRef]  

14. V. Akondi, C. Falldorf, S. Marcos, et al., Opt. Express 23, 25425 (2015). [CrossRef]  

15. K. Wang, Y. Li, Q. Kemao, et al., Opt. Express 27, 15100 (2019). [CrossRef]  

16. R. G. Paxman, T. J. Schulz, and J. R. Fienup, J. Opt. Soc. Am. A 9, 1072 (1992). [CrossRef]  

17. R. G. Paxman and J. Fienup, J. Opt. Soc. Am. A 5, 914 (1988). [CrossRef]  

18. R. Gerchberg and W. Saxton, Optik 35, 237 (1972).

19. R. A. Gonsalves, Opt. Eng. 21, 829 (1982). [CrossRef]  

20. H. Bao, C. Rao, Y. Zhang, et al., Opt. Lett. 34, 3484 (2009). [CrossRef]  

21. G. Wang and K. Wang, Optik 124, 6713 (2013). [CrossRef]  

22. Z. Li and X. Li, Opt. Express 26, 31675 (2018). [CrossRef]  

23. D. G. Voelz, Computational Fourier Optics: A MATLAB Tutorial (SPIE Press, 2011), p. 51.

Data availability

Data underlying the results presented in this paper are not publicly available but may be obtained from the authors upon reasonable request.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (5)
Fig. 1.
Fig. 1. Wrapped phases in the presence of different scales of noise and their corresponding HS spot arrays and PSFs. (a)–(c) Corresponding to the small scale, (d)–(f) corresponding to the middle scale, and (g)–(i) corresponding to the large scale.

View in Article | Download Full Size | PDF

Fig. 2.
Fig. 2. Principle of the PD-VHS algorithm.

View in Article | Download Full Size | PDF

Fig. 3.
Fig. 3. Sample of PSF to be filtered.

View in Article | Download Full Size | PDF

Fig. 4.
Fig. 4. Simulations: evaluation of quality-map guided, VHS, and PD-VHS algorithms for two scales of (a) $0.5 \lambda$ and (b) $1 \lambda$ aberrations with different scales of noise (SNR is shown in decibels).

View in Article | Download Full Size | PDF

Fig. 5.
Fig. 5. Example of PU results of quality-map guided, VHS, and PD-VHS algorithms. (a) Absolute phase (RMS = $1 \lambda$), (b) wrapped phase, (c)–(e) PU results of quality-map guided (RMSE = $4.91\lambda$), VHS (RMSE = $0.89\lambda$), and PD-VHS (RMSE = $0.07\lambda$) algorithms, (f)–(j) PSFs corresponding to phases (a)–(e).

View in Article | Download Full Size | PDF

Tables (1)

Tables Icon

Algorithm 1. GS algorithm

View Table

Equations (8)

Equations on this page are rendered with MathJax. Learn more.

φ ( x , y ) = ϕ ( x , y ) + 2 π k ( x , y ) + N ( x , y ) ,
ϕ ( x , y ) = ∑ i = 1 a i Z i ( x , y ) ,
P S F = | F ( P u p i l ( x , y ) ⋅ e ( j ⋅ ( ϕ ( x , y ) ) ) ) | 2 ,
  P S F = B u t t e r w o r t h _ f i l t e r [ P S F − T ] ( ( P S F < 0 ) = 0 )
T = α × m e a n 2 ( A + B + C + D ) + β × s t d 2 ( A + B + C + D ) .
A = D + G
[ a 1 a 2 ⋮ a v ] = [ Z 11 Z 12 ⋯ Z 1 v Z 21 Z 22 ⋯ Z 2 v ⋮ ⋮ ⋱ ⋮ Z 2 u , 1 Z 2 u , 2 ⋯ Z 2 u , v ] + [ g 1 g 2 ⋮ g 2 u ] ,
F W H M = λ ⋅ f D ⋅  pixel, 
Select as filters
Select Topics Cancel
Top
       
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.
Privacy | Terms of Use