# Robust processing of phase dislocations based on combined unwrapping and inpainting approaches

Haiting Xia, Silvio Montrésor, Rongxin Guo, Junchang Li, François

Olchewsky, Jean-Michel Desse, Pascal Picart

► To cite this version:

Haiting Xia, Silvio Montrésor, Rongxin Guo, Junchang Li, François Olchewsky, et al. Robust processing of phase dislocations based on combined unwrapping and inpainting approaches. Optics Letters, 2017, 42 (2), pp.322-325. 10.1364/OL.42.000322. hal-01445617

HAL Id: hal-01445617

https://hal.science/hal-01445617

Submitted on 2 Jul 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Robust processing of phase dislocations based
on combined unwrapping and inpainting
approaches

HAITING XIA,<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> SILVIO MONTRESOR,<fnref n="3" /> RONGXIN GUO,<fnref n="1" />,<fnref n="2" /> JUNCHANG LI,<fnref n="4" /> FRANÇOIS OLCHEWSKY,<fnref n="5" />
JEAN-MICHEL DESSE,<fnref n="5" /> AND PASCAL PICART<fnref n="3" />,\*

This Letter proposes a robust processing of phase disloca-
tions to recover continuous phase maps. The approach is
based on combined unwrapping and inpainting methods.
Phase dislocations are determined using an estimator based
on the second order phase gradient. The algorithm is
validated using a realistic simulation of phase dislocations,
and the phase restoration exhibits only weak errors. A
comparison with other inpainting algorithms is also pro-
vided, demonstrating the suitability of the approach. The
approach is applied to experimental data from off-axis
digital holographic interferometry. The phase dislocation
from phase data from a wake flow at Mach 0.73 are iden-
tified and processed. Excellent phase restoration can be
appreciated.

Phase retrieval and processing are unavoidable processes that
can be encountered in many applications for which the phase
is of primary interest, for example, synthetic aperture radar im-
aging [1], magnetic resonance imaging [2], surface shape mea-
surement [3], digital holography [4], and digital holographic
microscopy [5]. Phase processing requires recording a set of im-
ages and application of the arctangent operator to get a wrapped
(modulo 2π) phase. Then, unwrapping is a process of retrieving
the actual unwrapped phase from this wrapped phase. From a
practical point of view, phase unwrapping may be quite difficult
because of the presence of residues, noise, or other phase sin-
gularities. Noise can be processed with adapted spatial filtering
[6]. However, the existence of dislocations in the phase maps

will make phase unwrapping fail. Dislocations are related to
structural discontinuity in the phase fringe pattern. They may
be due to structural discontinuities or to singularities of the
measured field. For this kind of dislocations, the phase in the
dislocation region does not need any restoration. Unfortu-
nately, dislocations may be due to artifacts from noise reduc-
tion, or from under-sampling during recording (lack in spatial
or temporal resolution) and, generally, they occur in regions
where phase jumps are very close to one another. It follows that
the phase in the dislocation region needs to be restored. In prac-
tical phase maps, dislocations include not only noise, residues,
and discontinuities, but also continuous dislocations. In con-
ventional phase unwrapping algorithms, residues, quality maps,
and modulation are often used to deal with noise and disconti-
nuities [7]. In branch-cut algorithms, opposite residue charges
are connected by branch cuts that define the paths that cannot
be crossed in unwrapping. This approach can deal with some
dislocations such as residues, but cannot deal with regional
phase dislocations. A quality map is effective in dealing with the
noisy regions in the phase map, but its ability to process real
discontinuities is limited. Thus, conventional phase unwrap-
ping approaches can only deal with discontinuities and noise,
but not further with the continuous dislocations. In addition,
conventional unwrapping approaches do not restore any phase
dislocations such as branch cuts and low-quality regions. In this
Letter, a robust approach for efficiently processing phase dislo-
cations is proposed. The method is based on detecting dislo-
cations with a second-order phase gradient, masking regional
phase dislocations, least-squares unwrapping, and then inpaint-
ing the unwrapped phase.

Before processing, the dislocations need to be clearly iden-
tified in the phase map. Recently, Meng et al. [8] used modu-
lation to detect and mask out dislocations. Singularities were
processed by exemplar-based algorithm before unwrapping.
However, modulation calculation requires several phase fringe

<footnote n="1">¹Faculty of Civil Engineering and Mechanics, Kunming University of Science and Technology, Kunming 650500, China</footnote>

<footnote n="2">²Key Laboratory of Yunnan Higher Education Institutes for Mechanical Behavior and Microstructure Design of Advanced Materials, Kunming University of Science and Technology, Kunming 650500, China</footnote>

<footnote n="3">³Université du Maine, CNRS UMR 6613, LAUM, Avenue Olivier Messiaen, 72085 Le Mans Cedex 9, France</footnote>

<footnote n="4">⁴Faculty of Science, Kunming University of Science and Technology, Kunming 650500, China</footnote>

<footnote n="5">⁵ONERA, The French Aerospace Lab, 5Bd Paul Painlevé, 59000 Lille, France</footnote>

<footnote>*Corresponding author: pascal.picart@univ-lemans.fr</footnote>
patterns and cannot be performed with one single phase map. In addition, an exemplar-based algorithm is suitable for in-painting texture figures and small data-missing regions, but not for large dislocation area (several tens of pixels). Thus, this approach is not suitable for the phase map with large dislocations. In this Letter, segmentation of the phase map to identify and isolate large dislocation areas is based on the second-order gradient of the phase. Consider $\varphi_{ij}$ to represent the true phase and $\psi_{ij}$ to represent the wrapped phase at the grid point $(i, j)$ of a phase map, and note $W(\varphi_{ij}) = \psi_{ij}$ the wrapping operator [7]. The wrapped phase derivatives are defined as

$$ 
\begin{aligned}
\Delta_{ij}^x &= W(\psi_{(i+1)j} - \psi_{ij}), \\
\Delta_{ij}^y &= W(\psi_{i(j+1)} - \psi_{ij}),
\end{aligned} 
\qquad (1) 
$$

and the gradient of the true phase is expressed as

$$ G_{ij} = |\nabla \varphi_{ij}| = \sqrt{(\Delta_{ij}^x)^2 + (\Delta_{ij}^y)^2}. \qquad (2) $$

From Eqs. (1) and (2), the second-order wrapped phase derivatives are defined as

$$ 
\begin{aligned}
\Delta G_{ij}^x &= G_{(i+1)j} - G_{ip} \\
\Delta G_{ij}^y &= G_{i(j+1)} - G_{ij}.
\end{aligned} 
\qquad (3) 
$$

The second-order phase gradient can be computed with

$$ |\nabla G_{ij}| = \sqrt{(\Delta G_{ij}^x)^2 + (\Delta G_{ij}^y)^2}. \qquad (4) $$

For a continuous phase map, the values of the second-order phase gradient are weak. However, at the borders of the phase dislocations, its values are relatively large. Hence, it can be considered as a pertinent indicator to detect the borders of dislocations and, by choosing a threshold, a binary mask can be generated from its contour. Because the second-order phase gradient values of few points inside the dislocation area are minor, these points are set to zero in the mask. According to the residue theorem for phase unwrapping [7], the residues will cause path dependent results and make unwrapping fail. The sign of the residue is referred to as its polarity or charge. A mask can be used to connect the residues of opposite polarity or residues and the border of the phase map to avoid path-dependent results. In a masked region, the number of positive and negative residues should be equal, or the mask connects the residue and the border. A second-order phase gradient can detect residues, and automatic mask generation can be obtained using additional dilation. After applying the mask on the wrapped phase to mask out phase dislocations, the masked wrapped phase is unwrapped by a phase unwrapping algorithm based on least-squares and iteration. The algorithm is described in [9] and is referred to in this Letter as PULSI. This choice is related to a previously published paper [10] that demonstrates that its calibrated version is the most robust to noise.

An inpainting algorithm based on linear interpolation is used to restore the phase in the masked areas after unwrapping. The reason is that dislocations are mainly located in regions with close fringes, where the phase varies linearly. The two-dimensional interpolation of smooth images can be formulated as a variation problem of finding the function $\varphi(x, y)$ that minimizes the integral in the masked regions [11] [Eq. (5)], under the condition that $\varphi(x, y)$ reaches the given values along the contours. The calculation of the variation is used to find the conditions under which $J$ assumes a stationary value:

$$ J = \iint \left( \frac{\partial \varphi}{\partial x} \right)^2 + \left( \frac{\partial \varphi}{\partial y} \right)^2 dxdy. \qquad (5) $$

Therefore, the function $\varphi(x, y)$ that minimizes $J$ must be the solution of the partial differential equation that results from the appropriate Euler-Lagrange equation

$$ \frac{\partial^2 \varphi}{\partial x^2} + \frac{\partial^2 \varphi}{\partial y^2} = 0, \qquad (6) $$

under the condition that $\varphi(x, y)$ reaches the given values along the contours. Approximating Eq. (6), using finite differences for partial derivatives, and considering the unwrapped phase to be processed, implies that at any node in the grid, the inpainted value is given by the average of its four neighbors according to

$$ \varphi_{ij} = \frac{1}{4} [\varphi_{(i-1)j} + \varphi_{(i+1)j} + \varphi_{i(j-1)} + \varphi_{i(j+1)}]. \qquad (7) $$

If there are $n$ points in the masked region, one obtains $n$ linear equations. Thus, recovering missing data in the masked region is performed by solving the set of linear equations. In summary, the proposed approach to processing large dislocated phase areas is, first, detection and masking out dislocations with second-order phase gradients; second, unwrapping with PULSI; and, finally, inpainting by interpolation. This approach is referred to as PULSI+INTERP in the following.

In order to evaluate the robustness of the proposed approach, a numerical simulation was carried out to produce phase maps with realistic phase dislocations. For the simulation, a mathematical modeling with spatial variations is used to compute a phase (peaks Matlab function with 1024 × 1024 pixels). The number of modulo $2\pi$ fringes in the wrapped phase is controlled by the peak-to-valley. To generate dislocations, a strong moving average filter is applied to the sine and cosine of the phase. Then, after filtering, the phase is computed anew using the arctangent formula. In zones where the fringe density is high, compared to the size of the filtering kernel, phase dislocations are generated. To evaluate the performances of the proposed approach, the exact phase is computed by applying the same filtering to the original phase (that is, the continuous phase). With an initial phase at 1024 × 1024 pixels and a moving average filter of 25 × 25, one gets a simulated phase map of 975 × 975 with dislocations (pixels at the borders are removed). Figure 1(a) shows the original continuous phase, and Fig. 1(b)

<figure><img src="image_2.png" /><figcaption>Fig. 1. (a) Original phase, (b) wrapped phase with dislocations, (c) masked wrapped phase, and (d) unwrapped phase.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Inpainted phase with dislocation restoration from unwrapped phase (a) with PULSI+INTERP, (b) with PULSI+DIFF, (c) with PULSI+SRDCT, and (d) with PULSI+EB.</figcaption></figure>

shows the generated wrapped phase with dislocations. There are two large dislocation areas that can be observed in the upper lobe and the lower lobe of Fig. 1(b). As can be seen, the continuity of the fringes is broken, and unwrapping cannot be applied to these areas. When applying the proposed method, the phase dislocations can be selected and masked. When generating the mask of the dislocated phase maps, the threshold of the second-order phase gradient was set to 0.0045. This threshold was selected from the observation of the histogram of the second-order gradient. Figure 1(c) shows the masked areas. Figure 1(d) exhibits the unwrapped phase from PULSI, on which can be observed the masked zones. The last step is to apply inpainting to phase from Fig. 1(d). Figure 2(a) shows the final phase obtained from the full processing. At the naked eye, the result appears to be in very good agreement with the original phase from Fig. 1(a).

Note that there exist other approaches to inpaint image data. In order to demonstrate the suitability of the method proposed in this Letter for phase data, we compared the phase from Fig. 2(a) with other inpainting algorithms. Figure 2(b) shows the phase obtained after applying anisotropic diffusion (referred as PULSI+DIFF, with 100 iterations) [12]; Fig. 2(c) shows the phase obtained from sparse representation with over-complete discrete cosine transform dictionary [13] (referred as PULSI+SRDCT); and Fig. 2(d) shows the phase obtained with exemplar-based algorithm [14] (referred to PULSI+EB).

The quantitative appraisal of the phase quality obtained from these four inpainting approaches can be estimated by computing the phase error from the original unwrapped phase. Figure 3 exhibits the four error maps corresponding to the four restored phase maps from Fig. 2. As can be seen in Fig. 3, the error given by PULSI+INTERP is the lowest compared to the three other restorations. The quantitative evaluation of these error maps is provided in Table 1 which gives the global standard deviations $\sigma_{err}$, the standard deviations in the masked areas $\sigma_{errM}$, the global peak-to-valley error $PV_{err}$, and the peak-to-valley error in the masked areas $PV_{errM}$. From Table 1, it can be seen that standard deviation and peak-to-valley errors are the lowest with the proposed restoration approach. In addition, Table 1 gives the computation times ($t$) for the four-phase computation of Fig. 2. The PULSI+INTERP method exhibits the best performance.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Errors for restored phases (a) with PULSI+INTERP, (b) with PULSI+DIFF, (c) with PULSI+SRDCT, and (d) with PULSI+EB.</figcaption></figure>

<table><caption>**Table 1.** Statistics of Errors and Computation Time for the Four Different Inpainting Algorithms (PV is the Peak-to-Valley of the Unwrapped True Phase Map)</caption><thead><tr><th>Method</th><th>t(s)</th><th>σ<sub>err</sub><br>(rad)</th><th>σ<sub>errM</sub><br>(rad)</th><th>PV<sub>err</sub><br>(rad)</th><th>PV<sub>errM</sub><br>(rad)</th><th>PV<sub>err</sub>/PV<br>(%)</th></tr></thead><tbody><tr><td>PULSI+INTERP</td><td>27</td><td>0.061</td><td>0.386</td><td>2.01</td><td>1.85</td><td>3.4</td></tr><tr><td>PULSI+DIFF</td><td>165</td><td>0.209</td><td>1.000</td><td>5.19</td><td>4.31</td><td>8.9</td></tr><tr><td>PULSI+SRDCT</td><td>2665</td><td>0.146</td><td>1.236</td><td>6.15</td><td>6.15</td><td>10.5</td></tr><tr><td>PULSI+EB</td><td>169</td><td>0.127</td><td>1.073</td><td>5.66</td><td>5.66</td><td>9.7</td></tr></tbody></table>

The proposed approach was applied to experimental phase maps obtained from flow experiments. The setup is described in Fig. 4 which shows a digital off-axis holographic interferometer. The processing of the data is based on the two-dimensional Fourier transform to extract the optical phase [15]. The phenomenon of interest is a wake flow at Mach 0.73 propagating in the wind tunnel that is schematized by the glass windows in Fig. 4. A cylinder is inserted in the wind tunnel, and vortices are generated at upper side of the cylinder. The sensor is triggered, and the acousto-optic modulator (AOM) provides very short light pulses at 100 ns. However, at Mach 0.73, the exposure time provided by the pulse width may be not short enough to correctly record the high gradients of the refractive index. It follows that such high-speed transient variations generate phase dislocations in the phase maps due to the lack of temporal resolution. Figure 5 exhibits results obtained

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Experimental setup, AOM, and acousto-optic modulator; the wind tunnel is represented by the two glass windows.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Restoration of experimental phase maps by the proposed approach. (a) Raw wrapped phase (red circles indicate the dislocation areas), (b) masked wrapped phase, (c) unwrapped phase by PULSI, and (d) restored phase.</figcaption></figure>

from the processing of the phase. Figure 5(a) shows the wrapped phase from hologram processing. The red circles indicate areas in which phase dislocations can be observed. The masks generated by the proposed approach are shown in Fig. 5(b) with the uniform gray regions. The unwrapped phase from PULSI is given in Fig. 5(c). The final restored phase with INTERP is shown in Fig. 5(d). As can be seen, the continuity of the phase map in the discolorations areas is quite good. The effectiveness of the phase restoration can be appraised also from the cosine of the wrapped phase fringe pattern, as shown in Fig. 6(a).

Compared with the cosine fringes of the restored phase map, as shown in Fig. 6(b), one can appreciate that the restored fringes are in excellent agreement with those of the wrapped phase in the regions without dislocations. In the dislocation areas, the fringe continuity appears to be very well restored.

In conclusion, this Letter proposes a robust processing of phase dislocations to recover the continuous phase maps. The approach is based on combined unwrapping and inpainting methods. The processing requires detection and masking of the phase dislocations. The regions of interest are determined using an estimator based on the second-order phase gradient. The algorithm has been validated using a realistic simulation of the phase dislocations. Comparisons with three other inpainting algorithms demonstrated that the proposed method provides phase restoration with weak errors. In addition, the computation time was found to be the shortest of the set of tested algorithms. The approach was applied to the experimental data from off-axis digital holographic interferometry. Phase dislocation from the phase data from the wake flow at Mach 0.73 were identified and masked.

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** (a) Cosine fringe of original wrapped phase and (b) cosine fringe of restored phase.</figcaption></figure>

Then, unwrapping and inpainting yielded excellent phase estimation which was appreciated by analyzing the cosine of the phase fringe pattern. In the near future, the proposed method will be intensively applied to the phase restoration of the phase data from wake flow experiments up to Mach 2.0. Future work includes investigation of the robustness of the proposed approach to mix between strong speckle noise and phase dislocations.

**Funding.** National Natural Science Foundation of China (NSFC) (11362007, 11462009); Agence Nationale de la Recherche (ANR) (ANR-14-ASTR-0005-01).

## REFERENCES

1. R. Schlögel, C. Doubre, J. P. Malet, and F. Masson, *Geomorphology* **231**, 314 (2015).

2. P. Daga, T. Pendse, M. Modat, M. White, L. Mancini, G. P. Winston, A. W. McEvoy, J. Thornton, T. Yousry, I. Drobnjak, and J. S. Duncan, *Med. Image Anal.* **18**, 1132 (2014).

3. B. Tayebi, M. Jafarfard, F. Sharif, Y. Song, D. Har, and D. Kim, *Opt. Express* **23**, 11264 (2015).

4. P. Picart, *New Techniques in Digital Holography* (ISTE-Wiley, 2015).

5. P. Girshovitz and N. Shaked, *Opt. Express* **23**, 8773 (2015).

6. S. Montresor and P. Picart, *Opt. Express* **24**, 14322 (2016).

7. D. C. Ghiglia and M. D. Pritt, *Two-Dimensional Phase Unwrapping: Theory, Algorithms, and Software* (Wiley, 1998).

8. L. Meng, F. Sang, P. Yang, I. Wang, M. Komori, and A. Kubo, *Appl. Opt.* **51**, 2457 (2012).

9. H. T. Xia, R. X. Guo, Z. B. Fan, H. M. Cheng, and B. C. Yang, *Exp. Mech.* **52**, 439 (2012).

10. H. Xia, S. Montresor, R. Guo, J. C. Li, F. Yan, H. Cheng, and P. Picart, *Opt. Express* **24**, 28713 (2016).

11. S. Carlsson, *Signal Process.* **15**, 57 (1988).

12. P. Perona and J. Malik, *IEEE Trans. Pattern Anal. Mach. Intell.* **12**, 629 (1990).

13. M. Aharon, M. Elad, and A. Bruckstein, *IEEE Trans. Signal Process.* **54**, 4311 (2006).

14. A. Criminisi, P. Pérez, and K. Toyama, *IEEE Trans. Image Process.* **13**, 1200 (2004).

15. J. M. Desse, P. Picart, and P. Tankam, *Opt. Express* **16**, 5471 (2008).