Noised Phase Unwrapping Based on the Adaptive Window of Wigner Distribution

Abstract

A noised phase-unwrapping method is presented by using the Wigner distribution function to filter the phase noise and restore the gradient of the phase map. By using Poisson’s equation, the unwrapped phase map was obtained. Compared with the existing methods, the proposed method is theoretically simple, provides a more accurate representation, and can be implemented in light-field hardware devices, such as Shack-Hartmann sensors.

1. Introduction

Phase-unwrapping technology is key to restoring continuous phase distributions from measured phase maps [1]. It has been used in various fields, such as synthetic aperture radar (SAR) [2,3,4,5], structured light, atmospheric turbulence detection [6], the Gerchberg–Saxton algorithm [7], holographic interferometry [8,9,10], and surface measurement [11]. Owing to the periodic property of phase measurement, it only returns values in the range of (−π, π]. When the change in initial phase is greater than π, it will cause a discontinuity and jump. Phase-unwrapping technology is used to restore the original, smooth phase value with an interval of (−∞, ∞) from the discontinued wrapped phase map in the range of (−π, π].

To perform phase unwrapping, the branch-cut method [12,13], least-squares method [7,14], quality guide method [15,16], Fourier method [17], and many other methods [18,19,20,21,22] have been proposed. These procedures can be effectively performed when the phase map is correctly sampled, noise-free, and lacking abrupt phase changes due to the discontinuity in the original phase map. However, in experimental situations, the presence of noise changes the distribution of the wrapped phase map and leads to false discontinuities and jumps. These errors are transmitted to the surrounding areas and corrupt the entire unwrapping result.

To improve noise immunity, some prefiltering methods have been used to reduce the phase noise, such as the Kalman filter [23,24], adaptive phase filtering [25,26], and local frequency estimation [27,28,29]. Although these methods can effectively reduce the influence of noise, most of them are computer-based algorithms that are complicated, consume large amounts of computing time, and are difficult to connect with hardwired devices. It is necessary to introduce a new method that can not only filter the noise and perform phase unwrapping, but also connect with hardwired devices.

The development of light-field imaging and phase-retrieval technology has demonstrated the application prospects of wavefront detection and phase unwrapping using lenslet-based devices, such as the Shack–Hartmann (SH) sensor and light-field camera [30,31,32]. Those devices can be analyzed by phase-space optics theory, including the Wigner distribution function (WDF) [33]. The WDF was presented by Eugene Wigner in 1932 and then introduced into many other fields, such as signal analysis, light-field imaging, and image processing. It has the advantages of connecting the Fourier transform, describing both the position and direction information of optical rays, and connecting light-field imaging theory [33]. In addition to optical applications, it also has been used to estimate instantaneous frequencyand perform local image denoising and signal processing, which can be introduced into the optical field [29]. Thus, it is feasible and effective to perform phase unwrapping via WDF.

In this paper, a local phase-filtering and -unwrapping method based on gradient estimation by WDF is proposed. The phase map can be unwrapped by Poisson’s equation from the gradient of the phase map obtained from the WDF [34]. When the present method is applied into phase unwrapping based on an SH sensor, the phase map is unwrapped by Zernike modal recovery [35] from the gradient of the phase map obtained from the image captured by the SH sensor. The proposed procedure is path-independent and robust, can effectively reduce the effect of noise, and can perform phase unwrapping.

2. Theory of Unwrapping Based on WDF

The wrapped phase is defined by the relationship in Equation (1).

$$\phi \left(\overrightarrow{r}\right)=\psi \left(\overrightarrow{r}\right)-2\pi \cdot m\left(\overrightarrow{r}\right)$$

(1)

where $\overrightarrow{r}=\left(x,y\right)$ is the two-dimensional position vector, $\phi \left(\overrightarrow{r}\right)$ is the wrapped phase map, $\psi \left(\overrightarrow{r}\right)$ is the initial phase map, and $m\left(\overrightarrow{r}\right)$ is an array of integers to keep the range of the wrapped phase map within (−π, π]. After the wrapping process, the initial phase map with infinite range becomesa discontinuous, scoped wrapped phase, as shown in Figure 1.

Because it is difficult to perform information extraction and engineering applications on the wrapped phase map, owing to its discontinuity and limited boundary, an unwrapping method is necessary. Although there are huge differences between the original and wrapped phase map, they theoretically contain the same complex field from Equation (1), as shown in Equation (2).

$$P\left(\overrightarrow{r}\right)=\exp \left[i\cdot \phi \left(\overrightarrow{r}\right)\right]=\exp \left[i\cdot \psi \left(\overrightarrow{r}\right)\right]$$

(2)

where P is the complex field of the phase map. As WDF describes signals in both space and frequency simultaneously and independently [36,37], it is easier to fully analyze and process the frequency information of every single pixel. Furthermore, the connection between the WDF and optical devices also produces the processing method applied by light-field cameras and SH sensors, which are both widely used in atmospheric turbulence phase detection. Thus, the proposed method has the potential to perform phase unwrapping in both software and hardware. The WDF of a complex field of a phase map is established in Equation (3).

$$W_P(\overrightarrow{r},\overrightarrow{u})={\displaystyle \int {\Gamma }_P\left(\overrightarrow{r}+\frac{{\overrightarrow{r}}^{\prime }}{2},\overrightarrow{r}-\frac{{\overrightarrow{r}}^{\prime }}{2}\right)\exp \left(-i\overrightarrow{u}\cdot {\overrightarrow{r}}^{\prime }\right)d^2{\overrightarrow{r}}^{\prime }}$$

(3)

where $\overrightarrow{u}$ is the vector of the frequency domain coordinates, ${\Gamma }_P\left(\overrightarrow{r}+\frac{{\overrightarrow{r}}^{\prime }}{2},\overrightarrow{r}-\frac{{\overrightarrow{r}}^{\prime }}{2}\right)$ is the autocorrelation of the complex field of the phase map, which can be expressed as

$${\Gamma }_P\left(\overrightarrow{r}+\frac{{\overrightarrow{r}}^{\prime }}{2},\overrightarrow{r}-\frac{{\overrightarrow{r}}^{\prime }}{2}\right)=P\left(\overrightarrow{r}+\frac{{\overrightarrow{r}}^{\prime }}{2}\right)P^{\ast }\left(\overrightarrow{r}-\frac{{\overrightarrow{r}}^{\prime }}{2}\right)=\exp \left\{i\left[\psi \left(\overrightarrow{r}+\frac{{\overrightarrow{r}}^{\prime }}{2}\right)-\psi \left(\overrightarrow{r}-\frac{{\overrightarrow{r}}^{\prime }}{2}\right)\right]\right\}$$

(4)

By using the property of the Dirac function

$${\delta }^{\left(n\right)}\left(s\right)=\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }{\left(-ix\right)}^n}\exp \left(-isx\right)dx,\begin{array}{cc}& n=0,1,2\dots \end{array}$$

(5)

$${\displaystyle {\int }_{-\infty }^{\infty }f\left(x\right)}{\delta }^{\left(n\right)}\left(x\right)dx={\left(-1\right)}^n{f}^{\left(n\right)}\left(0\right),\begin{array}{cc}& n=0,1,2\dots \end{array}$$

(6)

where ${\delta }^{\left(n\right)}\left(s\right)$ is the nth-order derivatives of $\delta \left(s\right)$. By using the definition of the Wigner function and property of the Dirac function, the connection between the WDF and gradient of the phase map is established as

$$\nabla \psi \left(\overrightarrow{r}\right)=\frac{{\displaystyle \int W_P\left(\overrightarrow{r},\overrightarrow{u}\right)\overrightarrow{u}{d}^2\overrightarrow{u}}}{{\displaystyle \int W_P\left(\overrightarrow{r},\overrightarrow{u}\right)d^2\overrightarrow{u}}}$$

(7)

where $\nabla \psi \left(\overrightarrow{r}\right)$ are the gradients of the initial phase map. Equation (7) is the instantaneous frequency-estimation property of the WDF. Because the complex field of the initial phase map is the same as the wrapped phase map, the gradient of the initial phase map can be obtained by the wrapped phase map. If no singularities are present inside the initial phase map, the unwrapped phase map can be obtained by solving Poisson’s equation ${\nabla }^2\psi \left(\overrightarrow{r}\right)=0$ [38].

However, in engineering applications, the presence of noise affects the reconstruction of the phase map. To reduce the effect of noise, the filtering property of the WDF is used in this paper. The complex field of the noised phase map can be expressed as

$$P_n\left(\overrightarrow{r}\right)=P\left(\overrightarrow{r}\right)N\left(\overrightarrow{r}\right)=\exp \left[i\cdot \psi \left(\overrightarrow{r}\right)\right]\exp \left[i\cdot n\left(\overrightarrow{r}\right)\right]$$

(8)

where $P_n\left(\overrightarrow{r}\right)$ is the complex field of the noised phase map, $P\left(\overrightarrow{r}\right)$ is the complex field without noise and $N\left(\overrightarrow{r}\right)$ is the noise. The WDF of the noised phase map is defined by Equations (2) and (3).

$$W_{P_n}\left(\overrightarrow{r},\overrightarrow{u}\right)=W_P\left(\overrightarrow{r},\overrightarrow{u}\right){\ast }_{\overrightarrow{u}}{W}_N\left(\overrightarrow{r},\overrightarrow{u}\right)$$

(9)

where $W_{P_n}\left(\overrightarrow{r},\overrightarrow{u}\right)$, $W_P\left(\overrightarrow{r},\overrightarrow{u}\right)$, and $W_N\left(\overrightarrow{r},\overrightarrow{u}\right)$ are the WDFs of the complex field of the noised phase map, initial phase map, and noise, respectively. As $W_P\left(\overrightarrow{r},\overrightarrow{u}\right)$ has been discussed in Equation (3), $W_N\left(\overrightarrow{r},\overrightarrow{u}\right)$ should be analyzedforthe study of $W_{P_n}\left(\overrightarrow{r},\overrightarrow{u}\right)$. Because of its randomness, the noise cannot be directly acquired. By using Vina–Singin theorem, the Wigner–Ville spectrum [39] of the noise can be applied to describe the noise characteristic as

$$W_N(\overrightarrow{r},\overrightarrow{u})=={\displaystyle \int R_N\left[\overrightarrow{r}+\frac{{\overrightarrow{r}}^{\prime }}{2},\overrightarrow{r}-\frac{{\overrightarrow{r}}^{\prime }}{2}\right]\exp \left(-i\overrightarrow{u}\cdot {\overrightarrow{r}}^{\prime }\right)d^2{\overrightarrow{r}}^{\prime }}$$

(10)

where R_N is the autocorrelation function of the noise in a random process. In this paper, we assume that the noise is white noise with a range of [a, b]; thus, the distribution of R_N is [29]

$$R_N\left[\overrightarrow{r}+\frac{{\overrightarrow{r}}^{\prime }}{2},\overrightarrow{r}-\frac{{\overrightarrow{r}}^{\prime }}{2}\right]=\{\begin{array}{cc}1& {\overrightarrow{r}}^{\prime }=0\\ {\sin \mathrm{c}}^2\left(\frac{b-a}{2}\right)& {\overrightarrow{r}}^{\prime }\ne 0\end{array}$$

(11)

The Wigner spectrum of noise can be obtained from Equations (10) and (11) as

$$W_N(\overrightarrow{r},\overrightarrow{u})={\sin \mathrm{c}}^2\left(\frac{b-a}{2}\right)\delta \left(\overrightarrow{u}\right)+\left[1-{\sin \mathrm{c}}^2\left(\frac{b-a}{2}\right)\right]$$

(12)

Equation (12) contains two terms in the Wigner spectrum of the noise. The first term shows that there is a peak at zero frequency in the spectrum when the range of noise is less than 2π. The second term indicates that the power spectrum density is uniform in a nonzero frequency domain because Equation (12) is based on a random process that is statistically average. In instantaneous cases, the first term is a peak in the zero-frequency region while the second term is a random noise uniformly distributed across nonzero frequencies. From Equations (9) and (12), the WDF of the noised phase map can be described as

$$W_{P_n}\left(\overrightarrow{r},\overrightarrow{u}\right)={\sin \mathrm{c}}^2\left(\frac{b-a}{2}\right)W_P\left(\overrightarrow{r},\overrightarrow{u}\right)+\left[1-{\sin \mathrm{c}}^2\left(\frac{b-a}{2}\right)\right]$$

(13)

In Equation (13), the first term is the information of the initial phase map and the second term is the noise in the frequency domain, which should be filtered. Figure 2 shows the comparison of the WDF in the noised and initial phase map. Figure 2a,b show the initial and noised atmospheric turbulence phase maps. Figure 2c,d show the spatial frequency distribution of WDF at one pixel. From Figure 2c,d, it is easy to see that the frequency distribution of the noised phase map is constituted by the frequency distribution of the initial phase map and the noise.

It should be noted that the coordinate (184, 65) is a random point which is selected as an example to prove the local property of the WDF unwrapping method. The points can be replaced by any other points. From Equation (6), we can see that the gradients of Figure 2c,d are [0.05, 0.30] and [0.22, 2.01], respectively. Because the difference between these two gradients is large, the presence of noise in the frequency domain severely affects the gradient extraction. To reduce the effect of noise, filtering should be performed. In most cases, the frequency distribution of the initial phase map is centralized due to the smooth distribution of the initial phase, as shown in Figure 2c, while the frequency distribution of noise is uniform, as shown in Figure 2d. Thus, it is possible to use an adaptive window to separate the WDF of the initial phase map from that of the noised phase map via threshold segmentation. In this paper, the boundary of the windowis based on the intensity of the WDF. However, we cannot use the intensity of the initial WDF directly to establish the boundary of the window for the rapid change in the distribution, as shown in Figure 2d. To smooth the WDF, a high-frequency Gaussian filter is applied, namely the filtered WDF is

$$W_{p_2}\left(\overrightarrow{r},\overrightarrow{u}\right)={\displaystyle \int {\Gamma }_p\left(\overrightarrow{r},{\overrightarrow{r}}^{\prime }\right)}\exp \left(-\frac{{{\overrightarrow{r}}^{\prime }}^2}{{\sigma }^2}\right)\exp \left(-i\overrightarrow{u}\cdot {\overrightarrow{r}}^{\prime }\right)d^2{\overrightarrow{r}}^{\prime }$$

(14)

When the filtered WDF is added in Equation (7), we can obtain

$$\nabla \psi \left(\overrightarrow{r}\right)=\frac{{\displaystyle \int {\Gamma }_P\left(\overrightarrow{r},{\overrightarrow{r}}^{\prime }\right)\exp \left(-\frac{{{\overrightarrow{r}}^{\prime }}^2}{{\sigma }^2}\right)\exp \left(-i\overrightarrow{u}\cdot {\overrightarrow{r}}^{\prime }\right)\overrightarrow{u}{d}^2{\overrightarrow{r}}^{\prime }{d}^2\overrightarrow{u}}}{{\displaystyle \int {\Gamma }_P\left(\overrightarrow{r},{\overrightarrow{r}}^{\prime }\right)\exp \left(-\frac{{{\overrightarrow{r}}^{\prime }}^2}{{\sigma }^2}\right)\exp \left(-i\overrightarrow{u}\cdot {\overrightarrow{r}}^{\prime }\right)d^2{\overrightarrow{r}}^{\prime }{d}^2\overrightarrow{u}}}$$

(15)

We can see that the filtering process does not affect the accuracy of the phase-gradient extraction. Figure 3a,b are the filtered images of Figure 2c,d, respectively, and show that images become smoother after filtering.

From Equation (14), we find that the gradient estimation is [0.05, 0.29] in Figure 3a and [0.20, −2.02] in Figure 3b. Because the high-frequency filter smooths the spatial frequency distribution of both phase and noise, it is easier to establish the boundary of the adaptive window by reserving regions with high signal-to-noise ratios. However, when the ratio is small, it will affect the accuracy of boundary selection. Lower thresholds reserve the noise region, whereas larger thresholds lose the WDF of the initial phase map. According to Equation (13), because the noise is independent in both the spatial and frequency domains and the WDF of the initial phase map is correlative and continuous, we can superimpose the frequency distribution of the sampling point with adjacent pixels to enhance the weight of the phase map and reduce the weight of noise in the spatial frequency domain. However, this method inevitably increases the error of the gradient estimation by adding the frequency information of other pixels, as shown in Figure 4.

From Figure 4 and through comparison with Figure 3a,b, it can be seen that the weight of the noise is effectively reduced and the frequency distribution of the initial phase is retained because the gradient estimation of the phase map changed to [0.36, −0.70] from [0.20, −2.02]. Thus, the boundary of the window can be easily established.

Figure 5 is the filtered frequency distribution of Figure 4 obtained using image threshold segmentation, where (a) is the window depending on the intensity of the WDF and (b) is the processed distribution. It should be noted that only the central region is present in Figure 5 to clarify the window and processed distribution.

From Equation (6), we determine the gradient of Figure 5b as [0.04, −0.29], which is close to the original phase gradient at the sampling point. When the gradient of all pixels has been estimated, the unwrapped phase can be extracted by using Equation (7). Figure 6a shows the unwrapped phase map; Figure 6b shows the rewrapped phase map.

To validate the effectiveness of our method, the mean-squared error (MSE), which is defined by Equation (16), has been applied to describe the difference between the unwrapped and initial phase maps.

$$MSE=\frac{1}{M^2}{\displaystyle \sum {\displaystyle \sum Arg{\left\{\exp \left[i\left({\phi }_l\left(\overrightarrow{r}\right)-\phi \left(\overrightarrow{r}\right)\right)\right]\right\}}^2}}$$

(16)

where M is the size of the phase map and ${\phi }_l$ is the unwrapping phase. The MSE of the filtered phase map in Figure 5b is 0.033, while the MSE of the noised phase map in Figure 2b is 0.82. This shows that the noised, wrapped phase map can be effectively filtered and unwrapped using phase-gradient estimation based on WDF.

3. Phase Unwrapping by Imaging Processing

In this section, the performance of the proposed method is quantitatively and qualitatively evaluated with a spherical aberration phase map and a simulated atmospheric turbulence phase map. All phase maps contain noise with three different intensities. We compare our results with (1) the adaptive instantaneous frequency estimation proposed by Hussainand (2) the fast frequency-estimation approach proposed by Trouve.

3.1. Spherical Aberration Data

Figure 7 shows the wrapped spherical aberration phase map. To describe the effectiveness of our method, three different uniform noises have been added to Figure 8, as shown in the first column of Figure 8, where the next three columns show the results filtered by three different methods. The MSEs of the rewrapped phase maps are listed in

Table 1. Mean-squared error (MSE) of each filtering method in Figure 8.

Range of Noise	Noised Image	Instantaneous Frequency Estimation	Fast Frequency Estimation	Our Method
0.7π	0.404	0.0020	0.0020	0.0021
1.0π	0.824	0.0037	0.0064	0.0042
1.4π	1.606	0.025	0.07	0.014

.

3.2. Turbulence Phase Data

To confirm the effectiveness of the proposed method in a more complex application environment, such as atmospheric optics, an atmospheric turbulence phase map is generated. The turbulence phase map is shown in Figure 9. Three different uniform noises have been added to the first column of Figure 10, where the results filtered by three different methods are shown in the other three columns. The MSEs of the rewrapped phase maps are shown in

Table 2. MSE of each filtering method of Figure 10.

Range of Noise	Noised Image	Instantaneous Frequency Estimation	Fast Frequency Estimation	Our Method
0.7π	0.40	0.22	0.06	0.08
1.2π	1.18	0.21	0.15	0.08
1.4π	1.64	0.66	0.93	0.08

. From the figures and MSEs of the rewrapped phase maps, we can see that the proposed method can effectively reduce the phase noise while retaining details of the phase map. When the signal-to-noise ratio is lower, the proposed method can provide a better filtering effect than the other methods.

It should be noted that several pictures are tested in our work. Because all the examples give the same results, only the resultof one picture is shown in the paper.

4. Phase Unwrapping by Shack-Hartmann Sensor

In this section, to confirm the effectiveness of the proposed method by hardware system, the performance of the proposed method in an SH wavefront sensor is demonstrated using both simulated and experimental cases. The simulated SH sensor imaging algorithm was created by Sergio Bonaque-Gonzalez. The commercial SH sensor is Haso4 R-flex 50 made by Imagine Optics. Previous research has shown that the image captured by SH sensor can be considered as a smoothed WDF (SWDF) of the incident beam, as shown in Equation (17).

$$I\left(\overrightarrow{r}\right)=\frac{1}{\lambda f}{\displaystyle \int W_I\left({\overrightarrow{r}}^{\prime },{\overrightarrow{u}}^{\prime }\right)W_L\left({\overrightarrow{r}}^{\prime }-{\overrightarrow{r}}_0,{\overrightarrow{u}}^{\prime }-\overrightarrow{u}\right)d{\overrightarrow{r}}^{\prime }d{\overrightarrow{u}}^{\prime }}$$

(17)

where $I\left(\overrightarrow{r}\right)$ is the intensity distribution captured by a single lenslet on SH sensor, ${\overrightarrow{r}}_0$ is the center of a single lenslet, $\overrightarrow{u}$ is the spatial frequency, $W_I$ is the WDF of the incident beam, and $W_L$ is the WDF of the aperture of each lenslet. The SWDF is defined as the convolution of two WDFs over $\overrightarrow{r}$ and $\overrightarrow{u}$. Because the gradient estimation of $W_L$ is zero, we can obtain

$$\nabla \psi =\frac{{\displaystyle \int {\overrightarrow{u}}_1{W}_I\left({\overrightarrow{u}}_1\right)d{\overrightarrow{u}}_1}}{{\displaystyle \int W_I\left({\overrightarrow{u}}_1\right)}d{\overrightarrow{u}}_1}=\frac{{\displaystyle \int \overrightarrow{u}S\left[W_I,W_L\right]d\overrightarrow{u}}}{{\displaystyle \int S\left[W_I,W_L\right]d\overrightarrow{u}}}$$

(18)

where ${\overrightarrow{u}}_1={\overrightarrow{u}}^{\prime }$ and ${\overrightarrow{u}}_2=\overrightarrow{u}-{\overrightarrow{u}}^{\prime }$. Because smoothing $W_I$ does not affect the gradient estimation, the initial phase map can be restored by measuring $I\left(\overrightarrow{r}\right)$.

Then, the noise term is considered. As indicated in Equation (13), the noise is uniformly distributed in the spatial frequency domain. When the noise is smoothed by W_P, a smoothed, uniformly distributed noise is obtained, which is similar to background noise. Thus, we can reduce the phase noise by applying the same method that we used to reduce the background noise. Many methods have been used to reduce the background noise in wavefront reconstruction by SH sensor, some of which use threshold segmentation, which is the same as the method we have proposed in this paper [40]. This means that those methods can also be used to reduce the phase noise. By combining our method with existing denoising algorithms, we can improve its suitability for filtering phase noise.

In this experiment, the primary and noised 5th, 7th, 8–10th mixed, and 6th–10th mixed Zernike aberrations are introduced as the input wavefront with the spatial light modulator (SLM). Through phase-gradient estimation and Zernike modal recovering, the noised aberrations can be effectively unwrapped. In order to evaluate the performance of the proposed method both quantitatively and qualitatively, two traditional wavefront reconstruction algorithms are used for comparison with our method in both the simulation and experimental situations. The first algorithm is the full-spot method, which retains all captured data in the image; and the second is the edge-based threshold segmentation method, which filters the noise by referring to the noise distribution on the edge of $I\left(\overrightarrow{r}\right)$. To focus on both the image quality and optical properties of our method, both the MSE and Strehl ratio (STR) are considered.

4.1. Simulated Situation

Figure 11 shows the performance of noised phase unwrapping using traditional unwrapping methods and our method in the simulated situation. In a simulated, idealistic situation, in which there are no background noises, aberrations of the optical system, or other errors, only the phase noise is considered. We compare the results unwrapped by three different methods with the initial wavefront.

Table 3. MSE and Strehl ratio for each filtering method of Figure 10.

Zernike Coefficient	Initial Data		Without Denoising		Edge Threshold Denoising		Our Method
	MSE	STR	MSE	STR	MSE	STR	MSE	STR
5th Zernike Order	0.82	0.40	2.14	0.07	0.83	0.37	0.27	0.82
7th Zernike Order	0.82	0.40	3.22	0.16	0.60	0.48	0.30	0.73
8–10th Zernike Order	0.82	0.40	2.43	0.11	0.34	0.67	0.17	0.83
6–10th Zernike Order	0.82	0.40	2.61	0.04	0.29	0.74	0.16	0.82

shows the MSE of the noised phase map and unwrapping phase map. The results show that both methods can be used to filter phase noise, but our method has a better result.

4.2. Experimental Situation

Schematic for generating noised wrapped phase is shown in Figure 12.

The imagesfor different initial phasescaptured by SH sensor are given in Figure 13.

The SWDF of the noised phase map can be generated by using the images captured by the SH sensor. Then, we use the adaptive window to obtain the filtered gradient of the phase map and perform phase unwrapping. To show the differences between our methods with the traditional wavefront reconstruction method by SH sensor, the wrapped phases generated from these unwrapping phases are given in Figure 14. Through commercial SH sensors, the unwrapping phase maps can be experimentally restored. Because system errors are caused by background noise, aberrations of the optical system, errors of the SLM, and optical path errors of the incident beam, the detected phase map deviates from the input phase data. Thus, the unwrapping quality of noised Zernike aberrations is judged by the unwrapping results of the initial phase map, rather than by theoretical values.

Table 4. MSE and Strehl ratio for each filtering method of Figure 12.

Coefficient	Without Denoising		Edge Threshold Denoising		Our Method
	MSE	STR	MSE	STR	MSE	STR
5th Order	2.06	0.08	0.13	0.84	0.08	0.90
7th Order	2.59	0.10	0.10	0.88	0.10	0.89
8th–10th Order	2.94	0.13	0.22	0.77	0.19	0.80
6th–10th Order	2.83	0.04	0.24	0.74	0.18	0.79

shows the MSE of the noised phase map and unwrapping phase map. The results indicate that both filtering methods can be used in noise phase unwrapping.

5. Conclusions

A new approach based on WDF has been developed to filter noise and perform phase unwrapping on a noised phase map. Compared with existing methods, our method shows better filtering results on stronger noise for the local denoising property of WDF. To demonstrate the effectiveness of our method in image processing, we compare our results with the adaptive instantaneous frequency estimation proposed by Hussain and the fast frequency estimation approach proposed by Trouve. The results show that our method can provide a better result when the noise is stronger. We also demonstrate the effectiveness of the presented method in an SH sensor for the connection between WDF and light-field imaging. The results indicate that the distribution of phase noise is similar to that of the background noise and that our method can provide better results compared with traditional background denoising methods.