Deep Learning-Based Wrapped Phase Denoising Method for Application in Digital Holographic Speckle Pattern Interferometry

Abstract

This paper presents a new processing method for denoising interferograms obtained by digital holographic speckle pattern interferometry (DHSPI) to serve in the structural diagnosis of artworks. DHSPI is a non-destructive and non-contact imaging method that has been successfully applied to the structural diagnosis of artworks by detecting hidden subsurface defects and quantifying the deformation directly from the surface illuminated by coherent light. The spatial information of structural defects is mostly delivered as local distortions interrupting the smooth distribution of intensity during the phase-shifted formation of fringe patterns. Distortions in fringe patterns are recorded and observed from the estimated wrapped phase map, but the inevitable electronic speckle noise directly affects the quality of the image and consequently the assessment of defects. An effective method for denoising DHSPI wrapped phase based on deep learning is presented in this paper. Although a related method applied to interferometry for reducing Gaussian noise has been introduced, it is not suitable for application in DHSPI to reduce speckle noise. Thus, the paper proposes a new method to remove speckle noise in the wrapped phase. Simulated data and experimental captured data from samples prove that the proposed method can effectively reduce the speckle noise of the DHSPI wrapped phase to extract the desired information. The proposed method is helpful for accurately detecting defects in complex defect topography maps and may help to accelerate defect detection and characterization procedures.

1. Introduction

Protecting cultural relics and advocating for the protection of artworks with high artistic value is a mandatory aspect of cultural heritage management to secure the very existence of valuable cultural wealth for future generations. The timely inspection and repair of potential defects in cultural relics are of great significance for long-term preservation. In the last decades, many studies based on holographic interference have been proposed as a means for the implementation of new structural diagnosis tools, skills, and practices in heritage conservation [1,2,3,4,5,6,7,8,9,10,11,12]. Digital holographic speckle pattern interferometry (DHSPI) geometry has been proved to be an important technique applied in contactless art diagnostics non-destructive testing (NDT), since it provides high-resolution imaging of surface optical displacement, which can be used out of the laboratory for direct on-field non-destructive visualization and diagnosis of invisible structural alterations. Correspondingly, a fully portable DHSPI system providing high-density information and high fringe contrast data has been designed and developed for application in the remote examination of cultural heritage [13]. Methods of recovering data from speckle interferograms formed from another portable system based on digital speckle pattern interferometry were presented by Memmolo et al. [14].

Since artworks measured with DHSPI are usually complex, multilayered and inhomogeneous structures deteriorating with time and environmental conditions, starting from inside and affecting the precious surface [15,16,17,18,19], the resulting DHSPI wrapped phase is complicated, carrying complex fringe pattern changes and containing severe speckle noise. In many antique and heavily deteriorated artworks, the resulting interferograms may carry extremely discontinuous or even difficult to resolve fringe patterns. Noise can cause discontinuities in the fringes, making them appear to be interrupted. In this way, to properly process and evaluate the interferograms, tedious long-term fringe skeleton corrections and software for interactive interventions in the fringe pattern reconstruction process necessitating a priori knowledge are required [20,21]. Therefore, the denoising process is necessary in order to accurately observe and evaluate the state of the defects in artworks.

Many methods have been reported to denoise speckle noise from digital holography [22,23,24]. Phase-shifting technology provides the wrapped phase distributed in $\left(-\pi ,\pi \right]$, and a sine/cosine average filter (SCAF) [25] is proposed to denoise the numerator (sine fringes) and denominator (cosine fringes) of arctangent for processing noise and $2\pi$ discontinuities in the wrapped phase. Windowed Fourier transform (WFT) [26] is an effective method to denoise the wrapped phase, but it requires a long-processing time and adjusts the parameters to process data. The regularized quadratic cost function (RQCF) [27] is proposed by Jesús Villa et al. to denoise the wrapped phase, but this method requires the estimation of fringe orientation. There are still many efforts to automatically denoise the wrapped phase map while achieving high accuracy and fast processing. Recently, optical image denoising using deep learning has shown excellent processing speed and accuracy including fringe pattern denoising [28,29,30] and wrapped phase denoising [31]. Deep learning has demonstrated effective performance for reducing Gaussian noise in the wrapped phase [31], but the proposed method in [31] does not process the speckle noise of the wrapped phase for application in DHSPI.

In this paper, a deep learning method is further proposed for denoising DHSPI wrapped phase to help the structural diagnosis of artworks. In phase-shifting technology, the sine and cosine fringes are calculated by arctangent function to obtain the wrapped phase distributed in $\left(-\pi ,\pi \right]$. In order to address noise and $2\pi$ discontinuities of the wrapped phase, we present a deep learning method to denoise the sine and cosine fringes, and the filtered sine and cosine fringes are further calculated by arctangent to give an enhanced wrapped phase. Simulation results show that the proposed method can effectively reduce speckle noise and provide high-quality phase maps. The proposed method is verified by captured data from the portable DHSPI system. By comparing the high-frequency region of the phase map, the proposed method is shown to be superior to the post-processing of portable DHSPI systems and demonstrates a useful processing method for accurately assessing defects in artworks.

The rest of this paper is arranged as follows. Section 2 introduces the speckle fringe simulation method and phase retrieval method, which provide a theoretical basis for establishing a database. Section 3 presents the proposed method. Subsequently, Section 4 demonstrates the denoising performance from simulated data and real data. Section 5 concludes the paper.

2. Computer Simulation for Speckle Fringe and Phase Map

The schematic of the DHSPI system is shown in Figure 1, and the laser is divided into an object beam and a reference beam by a beam splitter. The object beam is driven by mirrors, and eventually illuminates the surface of the rough object. After the reference beam passes through a series of optical devices, the reference beam is combined with the object beam by a collimated beam combiner. The combined beam is recorded by a charge-coupled device (CCD) sensor. In the reference beam path, a mirror on a piezoelectric transducer is used for generating the optic path difference.

Assume that the object beam $U_O(u,v)$ and the reference beam $U_R(u,v)$ are described by

$$U_O(u,v)=a_O(u,v)exp\left[i{\varphi }_O(u,v)\right],$$

(1)

$$U_R(u,v)=a_R(u,v)exp\left[i{\varphi }_R(u,v)\right],$$

(2)

where $(u,v)$ are spatial coordinates, $a_O(u,v)$ and $a_R(u,v)$ are the corresponding amplitudes of beam, and ${\varphi }_O\left(u,v\right)$ and ${\varphi }_R(u,v)$ are the phase terms. Before the surface of the object is deformed, the resulting intensity $I_1(u,v)$ of the speckle pattern on the detector can be written as

$$I_1(u,v)={\left|U_O(u,v)+U_R(u,v)\right|}^2,$$

(3)

and we can further derive that

$$I_1(u,v)=a^2{}_O(u,v)+a^2{}_R(u,v)+2\sqrt{a^2{}_O(u,v)a^2{}_R(u,v)}\cos \left[{\phi }_1(u,v)\right],$$

(4)

where ${\phi }_1(u,v)={\varphi }_O(u,v)-{\varphi }_R(u,v)$. When the surface of the object is deformed, the resulting optical intensity $I_2(u,v)$ on the detector is defined as

$$I_2(u,v)=a^2{}_O(u,v)+a^2{}_R(u,v)+2\sqrt{a^2{}_O(u,v)a^2{}_R(u,v)}\cos \left[{\phi }_2(u,v)\right].$$

(5)

The phase difference between $I_1\left(u,v\right)$ and $I_2\left(u,v\right)$ can be expressed as

$$\mathsf{\Delta }\phi (u,v)={\phi }_1(u,v)-{\phi }_2(u,v),$$

(6)

and the subtraction operation between $I_1\left(u,v\right)$ and $I_2\left(u,v\right)$ can be expressed as a speckle fringe pattern,

$$I_d\left(u,v\right)=\left|I_1\left(u,v\right)-I_2\left(u,v\right)\right|.$$

(7)

The establishment of the database is necessary to train the model. In this work, the method of manipulating the spectrum of the light field by imitating the 4f system is used to generate a speckle pattern, and a low-pass filter $d(k_1,k_2)$ is simulated to be applied in object beams to imitate an aperture [32,33]. An object beam can be written as

$$U_O\left(u,v\right)=F^{-1}\left\{d(k_1,k_2)\times F\left\{a_O(u,v)exp[i{\varphi }_O(u,v)]\right\}\right\},$$

(8)

where $F$ represents the Fourier transform and F⁻¹ represents the inverse Fourier transform. For simulation, the phase ${\varphi }_{\mathrm{O}}(u,v)$ as a random variable is uniformly distributed in $\left(-\pi ,\pi \right]$, and the amplitude $a_O(u,v)$ is set to a constant of 1. The undeformed intensity $I_1\left(u,v\right)$ can be generated by Equation (3), where $U_R(u,v)=1$ for simulation. After the object deformation, the phase term ${\varphi }_{\mathrm{O}}(u,v)$ of the object beam is changed into ${\varphi }_{\mathrm{O}}(u,v)+\mathsf{\Delta }\phi (u,v)$. Similarly, the deformed intensity $I_2\left(u,v\right)$ can also be simulated by Equation (3). According to Equation (7), the speckle fringe pattern can be computed by the subtraction operation between $I_1\left(u,v\right)$ and $I_2\left(u,v\right)$.

The speckle fringe pattern with the different speckle sizes can be generated by adjusting the aperture size $d(k_1,k_2)$ according to Equation (8). An example with simulated data is shown in Figure 2, and the peaks function as the phase difference $\mathsf{\Delta }\varphi (u,v)$ is simulated. The speckle fringes (500 × 500 pixels) with the different aperture sizes generate different speckle sizes, where the smaller selected aperture produces a larger speckle size.

The wrapped phase can be estimated from the undeformed and deformed state of the intensity maps by a phase-shifting algorithm [34], and the wrapped phase is obtained by

$${\varphi }_N\left(u,v\right)=\mathrm{atan}2\left[S_N\left(u,v\right),C_N\left(u,v\right)\right],$$

(9)

where $S_N\left(u,v\right)$ are sine fringe patterns, and $C_N\left(u,v\right)$ are cosine fringe patterns. ${\varphi }_N\left(u,v\right)$ is the wrapped phase containing speckle noise, which is distributed in $\left(-\pi ,\pi \right]$. Figure 3 indicates the estimated wrapped phases under different aperture sizes. Here, the wrapped phase maps reveal the evident noise, which may produce errors in assessing defects of artworks. Therefore, effective filtering of the noisy wrapped phase map is important to help the structural diagnosis of artworks.

3. Deep Neural Network for DHSPI Wrapped Phase Denoising

3.1. Denoising Method

Convolutional neural networks (CNN) have been successful in processing various visual tasks and have demonstrated effective performance in handling Gaussian denoising [35]. In this work, the proposed network for DHSPI wrapped phase denoising is based on improved denoising convolutional neural networks (DnCNNs) [35]. The deeper neural network architecture for DHSPI wrapped phase denoising is illustrated in Figure 4.

In the proposed method, the noisy sine and cosine fringe patterns are denoising by CNN for solving the $2\pi$ discontinuities in the noisy wrapped phase. As can be seen from Figure 4, the noisy sine $S_N\left(u,v\right)$ and cosine $C_N\left(u,v\right)$ fringe patterns are input data. The input data are first computed by 35 convolutional layers, and the rectified linear unit (ReLU) [36] for non-linear mapping and batch normalization (BN) layers [37] for easier training are used after each convolutional layer. Each convolutional layer (from the 1st to 35th layer) has 32 filters of size 3 × 3 × 32. Supposing that the input data of each layer (from 1st to 35th) is ${\beta }_i$, the learned features $Z_i({\beta }_i)$ on each layer are expressed mathematically as [38]

$$Zi({\beta }_i)=ReLU\left[BN({\beta }_i*W_i)\right],\begin{array}{cccc}& & & i=1,\begin{array}{c}\end{array}2\begin{array}{c}\end{array}\dots \begin{array}{c}\end{array}35\end{array}$$

(10)

where $W_i$ represents the weights, and the * represents convolution operation.

The 36th layer is the convolutional layer, which outputs the denoised sine $S\left(u,v\right)$ and cosine $C\left(u,v\right)$ fringe patterns. Then, an enhanced wrapped phase map $\varphi \left(u,v\right)$ can be provided by arctangent function

$$\varphi \left(u,v\right)=\mathrm{atan}2\left[S\left(u,v\right),C\left(u,v\right)\right].$$

(11)

Network training is performed by finding the minimum loss function between the output data $O(k,r)$ and the ground truth $G(k,r)$. The adaptive moment estimation (Adam) optimizer [39] is chosen as the optimization method. A mini-batch size is expressed in $b$, and the used loss function $L$ is as follows

$$L=\frac{{\displaystyle \sum _{i=1}^b{\Vert O(k,r)-G(k,r)\Vert }^2}}{b}.$$

(12)

In this work, the learning rate is set to $1\times {10}^{-4}$. The weights are automatically updated after each iteration, and the training is completed after 100 epochs. Pytorch framework is used to perform network training, and the training process is processed on Ubuntu 14.04 with an Intel^® Core^TM i7-7820X [email protected] × 13 using the NVIDIA GTX 1080.

3.2. Data Set Preparation

To generate a data set for training the model, a simulation method of phase change $\mathsf{\Delta }\phi (u,v)$ is further presented here. The phase change $\mathsf{\Delta }\phi (u,v)$ is related to the out-of-plane displacement $H(u,v)$ as shown in Equation (13) [34], where $\lambda$ is wavelength (532 nm). First, the displacement $H(u,v)$ with size 80 × 80 pixels is simulated based on Gaussian function as shown in Equation (14). The arithmetic operations of addition and subtraction on 8 sets of Gaussian surfaces with the different mean value $\eta$ (from 1 to 80) and the different standard deviation $\epsilon$ (9, 13, and 18) are executed to obtain the displacement $H(u,v)$. Subsequently, phase changes with size 80 × 80 pixels are produced according to Equation (13), and 30,000 wrapped phases are generated using the simulation method mentioned in Section 2. Finally, the noisy sine and cosine fringe patterns and the corresponding noise-free data are obtained to train the model.

$$\mathsf{\Delta }\phi (u,v)=\frac{4\pi \times H(u,v)}{\lambda },$$

(13)

$$H(u,v)=600\times \frac{1}{2\pi {\epsilon }^2}\exp \left(-\frac{{(u-\eta )}^2+{(v-\eta )}^2}{2{\epsilon }^2}\right).$$

(14)

4. Evaluation Performance

4.1. Evaluation in Simulated Data

The denoising results of simulated data are demonstrated to evaluate the performance. In Figure 5a, the noisy data with the size of 500 × 500 pixels are simulated under the speckle radius of 1 pixel, and the enhanced wrapped phase using the proposed method manifests a significant improvement on the noisy data as shown in Figure 5b. The phase unwrapping of Figure 5b using the least-squares algorithm by the discrete cosine transform (DCT) [40] is implemented in Figure 5c, which has a low error distribution as shown in Figure 5d. The Mean Squared Error (MSE) comparison between the noisy data and denoised data is provided by

Table 1. Performance of the proposed method.

Figure	Noisy Data (MSE)	Denoising Data (MSE)
Figure 5a	3.9881	0.5823
Figure 6a	4.3608	0.8654

, and it can be seen that the noise is greatly reduced from the noisy wrapped phase.

By adjusting the size of the aperture, the simulated data in the speckle radius of 2 pixels are analyzed as shown in Figure 6a. The performance is demonstrated from the denoising result, and the denoising effect is illustrated by visualization in Figure 6b. In Figure 6c, the phase unwrapping is implemented by DCT for illustrating the denoising effect, and Figure 6d shows the low error distribution for illustrating the feasibility. In Table 1, the MSE comparison between the noisy data and denoising data indicates the good noise suppression ability of the proposed method. Besides, the calculation time for a wrapped phase with the size of 500 × 500 pixels is 3.45 s.

4.2. Evaluation in Captured Data

4.2.1. Data Acquisition

The precise monitoring and visualization of micro-deformation are important for understanding the structure of the defects [41]. To show the denoising effect to help the structural diagnosis, the captured data from the portable DHSPI system [13] are used to verify the performance. Figure 7 shows the schematic diagram of the measuring device [42], where IR lamps are used to provide the thermal excitation for the object. The object response after the thermal excitation is measured and recorded by the DHSPI system, and the personal computer (PC) can visualize results from the recorded data.

The object to be measured is marquetry as shown in Figure 8a. Figure 8b–f shows the estimated wrapped phase from the portable DHSPI system, and every wrapped phase represents a different state of the sample while the sample is cooling down after thermal excitation with IR lamps. The image in Figure 8b is the wrapped phase with a time interval of 12 s after the excitation. Subsequently, there are the wrapped phase maps with defect reaction at 53 s, 103 s, 151 s, and 218 s in a cooling down process after thermal excitation, respectively. However, we can observe obvious noise from these wrapped phase maps. In the high-frequency regions of the phase denoted by the red arrow, noise causes phase distortion. In such cases, it is difficult to directly estimate the depth information of the defect.

4.2.2. Denoising for Captured Data

Post-processing software of the portable DHSPI system can process the recorded results to improve the visualization for extracting useful information. The basic post-processing tools include denoising interferogram, phase unwrapping, displaying phase unwrapped image, and 3D surface deformation map, etc. [13]. The noisy data of Figure 8b–f are processed by the post-processing of the portable DHSPI system as shown in column (a) of Figure 9. Then, the proposed method is further used for denoising as shown in column (b) of Figure 9. From the visual inspection, the clean phase maps are obtained by the proposed method.

To clearly compare performance, in Figure 9 the enlarged phases in the red regions of column (a) and blue regions of column (b) are shown in columns (c) and (d), and the defect state can be observed more intuitively from the enhanced wrapped phase map using the proposed method. It can be seen that the proposed method can handle the high-frequency regions of the phase well, which shows it to be a useful post-processing method to estimate the depth of the defects in artworks. In the proposed method, the calculation time is 3.60 s for a captured wrapped phase map with size 1017 × 1030 pixels. Figure 10a,c shows the phase unwrapping results of the noisy data (Figure 8) and filtered data (Figure 9b) using DCT. The phase unwrapping effect is further illustrated by visualization of the contours in Figure 10b,d, and it further proves that the smooth phase can be obtained after denoising of the proposed method.

5. Conclusions

This paper presents a deep learning wrapped phase denoising method for application in digital holographic speckle pattern interferometry for the structural diagnosis of artworks. Based on simulation results, it is further verified that deep learning can be used to reduce speckle noise in the numerator (sine fringes) and denominator (cosine fringes) of the arctangent function, and an enhanced DHSPI wrapped phase using the arctangent function is provided from the denoising numerator and denominator. The captured data show that this method can effectively process the high-frequency regions of the wrapped phases compared with the DHSPI post-processing, thereby improving the accuracy of defect estimation in artworks. In future work, we will explore the defect prediction problem of artworks using the deep learning method. Meanwhile, we will consider applying deep learning to deformation detection of biological specimens in biomedicine [43,44,45,46]. Inhomogeneous surface properties significantly hinder reliable phase unwrapping and may produce a wrong quantitative evaluation for phase data. Deep learning may be a useful technique to solve these problems.