Skeleton extraction and phase interpolation for single ESPI fringe pattern based on the partial differential equations

Fang Zhang; Danyu Wang; Zhitao Xiao; Lei Geng; Jun Wu; Zhenbei Xu; Jiao Sun; Jinjiang Wang; Jiangtao Xi

doi:10.1364/OE.23.029625

1. Introduction

As a nondestructive technique, electronic speckle pattern interferometry (ESPI) can be used to measure the displacement of the rough surface of a specimen, which has been widely used for the whole-field measurement of displacements and their derivatives [1]. With ESPI, electronic speckle fringe patterns are obtained by processing the two speckles obtained before and after deformation of the tested surface. The phase of the fringe patterns carries the information related to the displacements, and hence it is extracted to yield the later. The commonly used approaches to obtain the phase are the phase-shifting method and the fringe centerline method. The traditional phase shifting method requires at least three speckle fringe patterns with certain phase differences. In the measurement of objects in fast motion or in a temporally unstable environment, it is difficult to take several phase-shifted interferograms in a short period of time. With the development of the phase shifting technology, a new spatial phase shifting method appears [2], which needs only one image before deformed and one image after deformed. However, it still requires more complicated optical system, thus not suitable for some applications. Therefore, it is always desired to extract the phase from a single fringe pattern, and the method based on fringe centerline extraction is an effective solution. The main idea is to detect the positions of the strongest and weakest stripes of the fringes (centerlines of the fringes), where the phase are integer multiples of $π$ , and fill the phase values in-between by means of interpolation. The traditional centerline extraction method consists of five steps, including preprocessing of the fringe pattern (e.g., filtering, enhancement), binarization, thinning to obtain the fringe centerline, assigning the skeleton order, and interpolation to get the whole-field phase values. Among these five steps, skeleton extraction and phase interpolation are particularly important and challenging.

There are mainly two types of methods for fringe skeleton extraction. One is fringe extreme tracking [3], and the other is thresholding binary-fringe and thinning [4]. However, the performance of these methods is dependent on the quality of the initial ESPI fringe image and the effectiveness of the pre-processing methods (such as filtering, enhancement etc.). Recently, some efforts have been made toward the computation of skeletons directly from a gray-scale image, which can avoid the disadvantage of the fringe information missing due to the inappropriate pre-processing [5,6 ].

Since the ultimate purpose of skeleton extraction is to get the whole-field phase values, phase interpolation plays a decisive role. The traditional methods include the nearest interpolation, bilinear interpolation and bi-cubic interpolation [7]. These approaches only consider the correlations of the neighbor points, thus ignoring the global features of image. Meanwhile, these methods usually suffer from the difficulty to identify the image edge and the phase jump, which will destroy the gradient information or lead to staircase characteristics. Besides these methods, Tang et al [8] proposed to train a Backpropagation (BP) Neural Network based on the phase values in the centerline for phase interpolation, while the learning speed of the neural network is slow. In addition, the interpolation capability of network depends on the structure and the training parameters (such as the number of neurons, the number of network lays and the training epochs), which are very difficult to be determined. At present, the neural network structure can only be selected based on experience due to the lack of unified and complete theoretical guidance. Wang et al [9] proposed the radial basis function (RBF) interpolation method to obtain unwrapped phase values based on a skeleton map. However, the extracted phase map is not very smooth. The further smooth operation maybe affects the precision of the extracted phase.

Over the past few years, partial differential equation (PDE) has demonstrated to be a powerful method and has been successfully applied in image filtering, enhancement, edge detection and segmentation. Since PDE based method is flexible, accurate and easy to operate, it is also applied in ESPI, such as filtering [10], binarization [11], skeleton extraction, etc. Jang and Hong [12] calculated a pseudo-distance-map using a nonlinear governing equation and extracted the skeletons from this pseudo-distance-map. Direkoglu et al [5] proposed to diffuse image in the dominance of direction normal to the feature boundaries and obtain the skeleton strength map by computing the mean curvature of level-sets. However, due to the existence of speckle noise, Direkoglu’s method is not suitable for ESPI fringe patterns. Tang et al [6] proposed a coupled PDE model for the gradient vector field (GVF) of ESPI fringe patterns, where the fringe skeleton is obtained via topological analysis of GVF after image diffusion. However, a threshold is required to obtain the skeleton of ESPI fringe in the process of image topological analysis. Since the threshold is not adaptive, inappropriate threshold will lead to discontinuity in the resulting skeleton lines.

In this paper, we propose to use PDEs for both ESPI skeleton extraction and phase interpolation. When extracting the fringe skeletons, we map the physical divergence to ESPI fringe images and take advantage of PDE for obtaining the fringe skeleton. The main process includes three steps, i.e. initialization of the GVF from the original image, GVF update by an anisotropic PDE, and skeleton extraction through divergence analysis on the updated GVF. The proposed method can be directly used on the gray-scale image without any preprocessing, and the fringe information can be well preserved. The proposed phase interpolation method also has two steps, i.e. order assignment for each skeleton, and the whole-field phase interpolation using the heat conduction equation. This method is based on heat diffusion theory, where the phase values of skeletons are regarded as the initial energy, which spreading to surroundings continuously with time (i.e., equation evolution). In the step of skeleton extraction for ESPI fringe, compared with a most existed skeleton extraction methods based on thresholding binary-fringe and thinning, our method does not need preprocessing, which can reserve the fringe information well. Also, compared with the method in Ref [6], which is also based on the gray fringe pattern, the proposed method is more robust because it does not need to set the threshold for the updated GVF. In the step of phase interpolation, the proposed phase interpolation method is faster and more accurate than BP Neural Network used in Ref [8].

The paper is organized as follows. Section II and Section III, the novel skeleton extraction method and the phase interpolation method are described in detail. Then, in Section IV, numerical solution method of PDEs is introduced. Next, skeleton extraction and phase interpolation results obtained by the proposed method are presented in comparison with other existing approaches in Section V. Finally, conclusion is given in Section VI.

2. Skeleton extraction for fringe patterns based on PDE

2.1 The initial gradient vector field

Let $I_{0} : R^{2} \to R$ represent a gray-level fringe pattern to be processed, and I ₀(x,y) is the gray-level of pixel (x,y). The initial GVF ${\vec{F}}_{0}$ of I ₀ can be obtained by

{\vec{F}}_{0} (x, y) = u_{0} (x, y) \vec{i} + v_{0} (x, y) \vec{j}

where

u_{0} (x, y) = \frac{\partial I_{0} (x, y)}{\partial x}

and

v_{0} (x, y) = \frac{\partial I_{0} (x, y)}{\partial y}

are the two components of

{\vec{F}}_{0} (x, y)

,

\vec{i}

and

\vec{j}

are the unit coordinate vectors [6]. For the discretized image, u ₀(x,y) and v ₀(x,y) for pixel (x,y) can be calculated by the difference of the gray values of the neighbouring pixels.

2.2 The diffusion of the gradient vector field

The GVF should be updated with the aim to smooth the noise. This function is similar to the cases of image smoothing [13]. In order to remove the noise and simultaneously keep the image edges, an anisotropic PDE is introduced to update u ₀(x,y) and v ₀(x,y) respectively.

A problem with the isotropic heat conduction equation is that it diffuses in all directions and hence may destroy the fringe edges of ESPI patterns in the process of denoising. To solve this problem, an inner orthogonal coordinate system is introduced based on the initial fringe pattern I ₀, which is characterized by the local gradient direction ( $ξ = \nabla I_{0} / | \nabla I_{0} |$ ) and the perpendicular-gradient direction ( $η = \nabla^{⊥} I_{0} / | \nabla I_{0} |$ ), respectively as shown in Fig. 1 .

Fig. 1 The inner orthogonal coordinate system based on the features of fringe patterns.

Download Full Size | PDF

Then, based on such coordinate system, u and v can be updated as the following anisotropic PDE respectively:

{\begin{matrix} \frac{\partial u}{\partial t} = u_{η η} + C u_{ξ ξ} u (x, y, 0) = u_{0} (x, y) \\ \frac{\partial v}{\partial t} = v_{η η} + C v_{ξ ξ} v (x, y, 0) = v_{0} (x, y) \end{matrix}

where u and v are the components of GVF to be updated with image evolution (i.e. iteration times increasing) and u ₀ and v ₀ are their initial values.

u_{ξ ξ}

,

u_{η η}

and

v_{ξ ξ}

,

v_{η η}

are the second order spatial derivatives of u and v along

ξ

and

η

, respectively. C is a positive constant to limit the diffusion degree along the gradient direction, which can better adjust the gradient vector field. The ESPI image is the statistical result of a large amount of speckles. Therefore the bigger C may result in excessive diffusion along

ξ

direction, leading to destroy of the fringe edges. To solve this problem, the parameter C is set to a small positive constant to reduce the diffusion along

ξ

direction. By restricting diffusion in the gradient direction of fringe patterns, good denoising can be achieved without destroying the fringe edges. In general, it should be less than 0.001.

In order to calculate $u_{ξ ξ}$ , $u_{η η}$ and $v_{ξ ξ}$ , $v_{η η}$ in the inner orthogonal coordinate system, the orientation angle $θ$ of fringe should be introduced, which is defined as the angle between fringe orientation and the x axis as shown in Fig. 1. $θ (x, y)$ for each pixel (x,y) is estimated within a neighborhood of a small window $Α$ with its center on the current point (x,y) [14]. The window size is chosen based on the quality of the fringe patterns. That is due to the orientation obtained by the gradient method is easily affected by high speckle noise. As the window size increasing, the computation error decreases gradually. However, the large window size will cause the increase of computation time. Generally, the appropriate window size is between15 × 15 and 55 × 55 pixels [15].

θ (x, y) = \frac{1}{2} arc \tan [\frac{\sum_{(i, j) \in Α} 2 u_{0} (i, j) v_{0} (i, j)}{\sum_{(i, j) \in Α} (u_{0}^{2} (i, j) - v_{0}^{2} (i, j))}]

Hence

u_{ξ ξ}

and

u_{η η}

can be calculated by

{\begin{matrix} u_{ξ ξ} = u_{y y} \cos^{2} θ + u_{x x} \sin^{2} θ - 2 u_{x y} \sin θ \cos θ \\ u_{η η} = u_{x x} \cos^{2} θ + u_{y y} \sin^{2} θ + 2 u_{x y} \sin θ \cos θ \end{matrix}

v_{ξ ξ}

and

v_{η η}

can be calculated in a similar way. With u and v updated by (2), GVF

\vec{F}

can be obtained as follows:

\vec{F} (x, y) = u (x, y) \vec{i} + v (x, y) \vec{j}

2.3 Skeleton extraction based on the divergence analysis

Divergence measures the magnitude of the source or sink of a vector field at a given point (x, y). The divergence of a continuously differentiable vector field $\vec{A} (x, y) = P (x, y) \vec{i} + Q (x, y) \vec{j}$ is given by:

\nabla \cdot \vec{A} (x, y) = \frac{\partial P (x, y)}{\partial x} + \frac{\partial Q (x, y)}{\partial y}

When the divergence of a point is greater than zero, i.e. $\nabla \cdot \vec{A} (x, y) > 0$ , the point has a positive source which emits flux in the vector field. When $\nabla \cdot \vec{A} (x, y) < 0$ , this point has a negative source which absorbs the flux. While the divergence is equal to zero means that there is no source at the point (x, y).

Let us consider the divergence of the GVF $\vec{F}$ as follows:

φ (x, y) = \nabla \cdot \vec{F} (x, y) = \frac{\partial u (x, y)}{\partial x} + \frac{\partial v (x, y)}{\partial y}

The above can be used to determine the characteristic of the point (x,y), as a source point or a sink point. When the divergence of the point (x,y) is greater than zero, (x,y) is a source point. If the divergence is less than zero, (x,y) is a sink point. By means of thresholding to

φ (x, y)

, two binary images can be obtained:

φ_{D} (x, y) = {\begin{matrix} 1 φ (x, y) > 0 \\ 0 φ (x, y) < 0 \end{matrix}, φ_{B} (x, y) = {\begin{matrix} 0 φ (x, y) > 0 \\ 1 φ (x, y) < 0 \end{matrix}

where

φ_{D} (x, y)

and

φ_{B} (x, y)

corresponding to the dark and bright areas respectively. The single pixel-wide skeleton for the dark fringe can be obtained through thinning

φ_{D} (x, y)

(e.g. morphologic thinning [16]). Similarly, the skeleton for the bright fringe can be extracted by thinning the binary image

φ_{B} (x, y)

. Such an idea can be implemented as follows:

Step 1: Calculate the initial GVF of the given fringe pattern ${\vec{F}}_{0}$ using Eq. (1);
Step 2: Calculate the orientation angle $θ$ for the initial ESPI fringe pattern using Eq. (3);
Step 3: Implement the anisotropic diffusion PDE in Eq. (2) on the two components u ₀ and v ₀, and obtain the updated GVF $\vec{F}$ ;
Step 4: Calculate the divergence of $\vec{F}$ using Eq. (7);
Step 5: Get the fringe skeletons according to the divergence property of $\vec{F}$ by means of thresholding via Eq. (8);
Step 6: Thinning using mathematical morphology thinning method.

3. Phase interpolation for fringe patterns based on PDE

3.1 Assignment of the order number of skeletons

With the skeleton images extracted, the next task is to allocate a sequential order number to each of the skeletons. The number is assigned in such a way that, if a particular skeleton is numbered as n, the adjacent one in the peak direction will be n + 1, and the adjacent one in the valley direction will be n-1. The phase on the fringe centerlines is the product of π by the skeleton orders. Let I _s denote the binary skeleton image of the fringe pattern, and the assigned skeleton image is denoted by I _a. (The sub-indexes “s” and “a” are the abbreviations of “skeleton” and “assigned”, respectively.)

3.2 Phase interpolation

We propose to interpolate the whole-field phase values based on the heat conduction theory. The phase image domain is regarded as a heat-conducting plane, with the phase values in the image being the heat (or temperature) distribution on the plane. The phase values on the fringe centerlines are known and can be seen the source of heat, which will diffuse over the surrounding area, causing the phase change of the whole field.

The heat conduction equation and the initial condition can be expressed as follows:

\partial_{t} I_{p} = \nabla^{2} I_{p}, I_{p} (x, y, 0) = I_{a} (x, y)

where I _a(x,y) is the initial phase of the point (x,y), I _p(x,y,t) is the interpolated phase of the point (x,y) at time t (“p” is the abbreviation of “phase”), whose initial value I _p(x,y,0) is I _a(x,y). In the heat diffusion process, energy diffuses from fringe skeletons to the whole field. As the energy of the skeletons will decrease with the diffusion, the phase values on the skeletons will decrease and deviate from the initial real phase values. Therefore, phase compensation should be implemented to make the phase values of the skeleton constant during the diffusion, that is:

I_{p} (x_{s}, y_{s}, t_{n}) = I_{a} (x_{s}, y_{s}), (x_{s}, y_{s}) \in S

where S is the set of the points on the skeleton. According to the above modified energy conditions, the interpolated whole-field phase of the fringe pattern can be obtained based on Eq. (9).

4. Numerical solution of PDEs

Equations (2) and (9) are continuous and should be discretized in order to have numerical solutions for updating the variables. To this end, the grayscale image to be processed is represented by a matrix of M × N intensity values, i.e, u(i, j) = u_i,j for 1<i<M, 1<j<N. With the PDEs, the image is updated at time instances $t_{n} = n Δ t$ , we denote u(i, j, t_n) by $u_{i, j}^{n}$ , where the subscripts i, j give the position in the two-dimensional grid, and the superscript n is the time index.

The time derivative u_t at (i, j, t_n) is approximated by the forward difference

{(u_{t})}_{i, j}^{n} = \frac{u_{i, j}^{n + 1} - u_{i, j}^{n}}{Δ t}

where

Δ t

is step-size.

The spatial derivatives can also be estimated as follows:

{(u_{x x})}_{i, j}^{n} = u_{i + 1, j}^{n} - 2 u_{i, j}^{n} + u_{i - 1, j}^{n}

{(u_{y y})}_{i, j}^{n} = u_{i, j + 1}^{n} - 2 u_{i, j}^{n} + u_{i, j - 1}^{n}

{(u_{x y})}_{i, j}^{n} = (u_{i + 1, j + 1}^{n} - u_{i - 1, j + 1}^{n} - u_{i + 1, j - 1}^{n} + u_{i - 1, j - 1}^{n}) / 4

In a similar way we can also determine

{(v_{x x})}_{i, j}^{n}

,

{(v_{y y})}_{i, j}^{n}

and

{(v_{x y})}_{i, j}^{n}

.

Based on the above Eq. (2) can be rewritten in the following discrete form, which can be used to update gradient vector field with every time increment:

{\begin{cases} u_{i, j}^{n + 1} = u_{i, j}^{n} + Δ t [{(u_{η η})}_{i, j}^{n} + C {(u_{ξ ξ})}_{i, j}^{n}] \\ v_{i, j}^{n + 1} = v_{i, j}^{n} + Δ t [{(v_{η η})}_{i, j}^{n} + C {(v_{ξ ξ})}_{i, j}^{n}] \end{cases}

where

{\begin{cases} (_{u_{η η}}^{)} = (_{u_{x x}}^{)} \cos^{2} (θ_{i, j}) + 2 (_{u_{x y}}^{)} \cos (θ_{i, j}) \sin (θ_{i, j}) + (_{u_{y y}}^{)} \sin^{2} (θ_{i, j}) \\ (_{u_{ξ ξ}}^{)} = (_{u_{y y}}^{)} \cos^{2} (θ_{i, j}) - 2 (_{u_{x y}}^{)} \cos (θ_{i, j}) \sin (θ_{i, j}) + (_{u_{x x}}^{)} \sin^{2} (θ_{i, j}) \\ (_{v_{η η}}^{)} = (_{v_{x x}}^{)} \cos^{2} (θ_{i, j}) + 2 (_{v_{x y}}^{)} \cos (θ_{i, j}) \sin (θ_{i, j}) + (_{v_{y y}}^{)} \sin^{2} (θ_{i, j}) \\ (_{v_{ξ ξ}}^{)} = (_{v_{y y}}^{)} \cos^{2} (θ_{i, j}) - 2 (_{v_{x y}}^{)} \cos (θ_{i, j}) \sin (θ_{i, j}) + (_{v_{x x}}^{)} \sin^{2} (θ_{i, j}) \end{cases}

and the boundary conditions of GVF diffusion are defined as follows:

u_{i, 0}^{n} = u_{i, 1}^{n}, u_{i, N + 1}^{n} = u_{i, N}^{n}, v_{i, 0}^{n} = v_{i, 1}^{n}, v_{i, N + 1}^{n} = v_{i, N}^{n}, i = 1, 2, \dots, M

u_{0, j}^{n} = u_{1, j}^{n}, u_{M + 1, j}^{n} = u_{M, j}^{n}, v_{0, j}^{n} = v_{1, j}^{n}, v_{M + 1, j}^{n} = v_{M, j}^{n}, j = 1, 2, \dots, N

Similarly, Eq. (9) can be discretized to yield numerical solutions to the whole-field phase interpolation by the following:

{(I_{p})}_{i, j}^{n + 1} = {(I_{p})}_{i, j}^{n} + Δ t ({(I_{p})}_{i + 1, j}^{n} + {(I_{p})}_{i - 1, j}^{n} + {(I_{p})}_{i, j + 1}^{n} + {(I_{p})}_{i, j - 1}^{n} - 4 {(I_{p})}_{i, j}^{n})

and the boundary conditions are defined by:

{(I_{p})}_{i, 0}^{n} = {(I_{p})}_{i, 1}^{n}, {(I_{p})}_{i, N + 1}^{n} = {(I_{p})}_{i, N}^{n}, i = 1, 2, \dots, M

{(I_{p})}_{0, j}^{n} = {(I_{p})}_{1, j}^{n}, {(I_{p})}_{M + 1, j}^{n} = {(I_{p})}_{M, j}^{n}, j = 1, 2, \dots, N

5. Experiments and results

In order to test the performance of the proposed approaches, experiments are carried out using a personal computer with Intel-core-i3 2.53 GHz, 2.00 GB of RAM and MATLAB R2012b.

For better illustration of the proposed skeleton extraction method, Fig. 2 shows the comparison results of our method and the method proposed in Ref [6]. In Figs. 2(a)-2(c) come from Ref [6]. Figure 2(a) is an initial experimentally obtained fringe pattern, which depicts the out-of-plane displacement of a circular plate. The plate is rigidly clamped at its boundary and is subjected to a central load. Figure 2(b) shows the black fringe skeleton based on the topological analysis of the GVF in Ref [6]; Fig. 2(c) gives the superimposition of the white fringe skeleton onto the initial ESPI fringe image; Fig. 2(d) and 2(e) are the white fringe skeleton and the black fringe skeleton obtained by our method with $Δ t = 0.2$ , $n = 300$ , $C = 0.0001$ and the window size to calculate orientation angle is 43 × 43; Fig. 2(f) gives the superimposition of Figs. 2(d) and 2(e) onto the initial ESPI fringe image. The results of the experiments show that the skeleton lines extracted by our method are consistent with the initial image and do not suffer from fracture phenomenon. However, as shown in Fig. 2(c) the white skeleton line is obviously fractured, which due to thresholding during skeleton decision in the method proposed by Ref [6]. The proposed method is better than the method in Ref [6]. in skeleton extraction.

Fig. 2 Experimentally-obtained ESPI fringe pattern (which comes from Ref [6].) and its skeletons. (a) Initial image; (b) Black fringe skeletons of (a) in Ref [6]; (c) Superimposition of white fringe skeletons of (a) on to (a) in Ref [6]; (d) White fringe skeletons of (a) by our method; (e) Black fringe skeletons of (a) by our method; (f) Superimposition of the white and black fringe skeletons onto (a).

Download Full Size | PDF

The next step is to assign the skeleton orders. Firstly, we superimpose the white fringe skeleton [Fig. 2(d)] onto the black fringe skeleton [Fig. 2(e)]. Then the skeletons are numbered from the outermost one to the inner ones on sequential basis (that is, 1, 2,…). We then multiply the skeleton order number by π, resulting in the initial heat distribution for energy diffusion as shown in Fig. 3(a) . With the proposed interpolation method, the phase image of Fig. 2(a) is obtained in Fig. 3(b), where obtained by the heat conduction equation [Eq. (9)] was updated 700 times with step size $Δ t = 0.25$ . In comparison, the phase image computed by BP Neural Networks (introduced by Ref [8].) is shown in Fig. 3(c). According to Ref [8], the locations of the pixels on the skeletons and the corresponding phase values in Fig. 3(a) are used as the training data sets in the BP Neural Networks for interpolation. In this test, the two hidden layers have 50 and 2 neurons, respectively. The network is trained for up to 1000 epochs to an error goal of 0.00001. In order to show the difference of the phase map obtained by these two methods, the three-dimensional phase graphs corresponding to Figs. 3(b) and 3(c) are presented in Figs. 3(d) and 3(e), respectively. It can be seen that the phase map obtained by our method is smoother than that computed by BP Neural Networks, especially at the edge of the whole field (see in Figs. 3(d) and 3(e)).

Fig. 3 Order assignment and phase interpolation for an experimentally-obtained fringe pattern. (a) Order assigned image; (b) Gray-phase image of (a) by our method; (c) Gray-phase image of (a) by BP Neural Networks; (d) and (e) The three-dimensional phase graph of (b) and (c).

Download Full Size | PDF

In order to quantitatively analyze advantages of our proposed method, we compared the accuracy and running time of the two interpolation methods. Figure 4(a) shows a simulated ESPI fringe pattern; Fig. 4(b) is the image of order-assigned skeletons obtained by the proposed skeleton extraction with $Δ t = 0.2$ , $n = 105$ , $C = 0.0001$ and the window size of orientation angle is 35 × 35; Fig. 4(c) gives the three-dimensional real phase graph used for simulation; Fig. 4(d) displays the pseudo-color image of the real phase; Fig. 4(e) is the pseudo-color phase image obtained by our interpolation method with $Δ t = 0.25$ , $n = 1000$ ; Figs. 4(f)-4(h) show the pseudo-color phase images computed by BP Neural Network interpolation method three times. Although the three results are obtained with the same network parameters, i.e. the hidden layers have 50 neurons, and the network is trained for up to 1000 epochs to an error goal of 0.00001, the phase images are different. The reason is that the network weights should be initialized with random values in the process of neural network training, which leads to different results under the same network parameters. In other words, the phase is not unique when using BP Neural Network interpolation even if the network parameters are fixed. While there is no such a problem in our PDE based interpolation method.

Fig. 4 Order assignment and phase interpolation for a computer-simulated fringe pattern. (a) A simulated original speckle fringe image; (b) Order assigned image for the skeletons of (a); (c) Three-dimensional graph of real phase; (d) Pseudo-color image of real phase; (e) Pseudo-color phase image by our method; (f), (g) and (h) are the pseudo-color phase image by BP Neural Network interpolation method (which come from three network training with the same parameters, respectively).

Download Full Size | PDF

In addition to uniqueness of the interpolation results, there are other two advantages of the proposed phase interpolation method. Firstly, the proposed method only needs 3.72s to interpolate the phase value [Fig. 4(e)], while the BP Neural Network requires the training time of 7s, 9s and 14s respectively corresponding to Figs. 4(f)-4(h). The reason of the different training time is the difference in the number of training iterations, also caused by random initiation. Secondly, the phase images obtained by BP Neural Network interpolation method are not smooth, especially a large area of phase fluctuation located at image boundary (see Figs. 4(f)-4(h)). While the phase image shown in Fig. 4(e) is smooth and very similar to the real phase image. Compared with the real phase value, the average phase errors associated with the proposed method and BP Neural Network interpolation are 0.01%, 3.28%, 4.37% and 4.41%, respectively. The reason for large error of BP Neural Network interpolation is that all training data are from the skeleton points, which can yield good result near the skeletons, but bad near the boundary due to the lack of sample information. The proposed method is based on the principle of heat conduction, which can extend the skeleton phase information to the boundary of image.

Besides, we also compared the phase extraction method reported in Ref [9]. and the proposed method. In Ref [9], the fringe extreme tracking method is adopted to extract the fringe skeletons. Based on the fringe skeletons, Ref [9]. gives a good interpolation. We give the comparison results based on an initial ESPI fringe pattern from Ref [9], which are shown in Fig. 5 . Figure 5(a) is the initial ESPI fringe pattern from Ref [9]. Figures 5(b) and 5(c) are the two skeleton maps from Ref [9]. and by the proposed skeleton extraction method. The interpolation results based on Fig. 5(b) by Ref [9]. and the proposed interpolation method are given in Figs. 5(d) and 5(e) respectively. The last image show the interpolation result based on Fig. 5(c) by the proposed interpolation method. Figure 5(c) is obtained with $Δ t = 0.2$ , $n = 650$ , $C = 0.0001$ and the window size of orientation angle is 43 × 43; Figs. 5(e) and 5(f) are achieved with $Δ t = 0.25$ , $n = 6000$ and $Δ t = 0.25$ , $n = 5000$ , respectively. It can be seen that based on the same skeleton map, our interpolation result [Fig. 5(e)] is smoother than the phase map [Fig. 5(d)] in Ref [9]. Synthesizing the proposed skeleton extraction method and the phase interpolation method based on PDE, the better phase map [Fig. 5(f)] can be obtained.

Fig. 5 Phase extraction results. (a) Initial ESPI fringe pattern (comes from Ref [9].); (b) Skeleton map of bright fringes from Ref [9]; (c) Skeleton map of bright fringes by our method; (d) Gray image of evaluated phase from Ref [9]. based on (b); (e) Gray image of evaluated phase obtained by our method based on (b); (f) Gray image of evaluated phase obtained by our method based on (c).

Download Full Size | PDF

In order to further demonstrate the performance of proposed approaches for skeleton extraction and phase interpolation, we computed the phase value of the dense and sparse computer-simulated fringe pattern, respectively. The results are shown in Fig. 6 . Figures 6(a-1) and 6(a-2) give a dense and a sparse simulated initial fringe pattern, respectively; Figs. 6(b-1) and 6(b-2) are the order-assigned images corresponding to the skeletons of Figs. 6(a-1) and 6(a-2) where $Δ t = 0.2$ , $C = 0.0001$ , $n = 40$ and $n = 200$ , the window size are 21 × 21 and 27 × 27, respectively; Fig. 6(c-1) is the three-dimensional phase graph of Fig. 6(a-1) by our method with $Δ t = 0.25$ , $n = 1000$ . It takes 2.86s to interpolate the phase value, and the error of phase value is 0.15% compared with the real phase value; Fig. 6(c-2) gives the three-dimensional phase graph of Fig. 6(a-2) by our method with $Δ t = 0.25$ , $n = 2000$ . The running time is 5.55s for the process of phase interpolation, and the error of phase value is 0.38% compared with the real phase value. From Fig. 6, it can be seen that our proposed skeleton extraction method can extract continuous and uncrossed skeletons even if the initial fringe pattern is very dense (see Fig. 6(b-1)). Also, the proposed interpolation method can give the smooth and exact whole-field phase even though the skeletons are sparse and have a little initial information (see Fig. 6(c-2)).

Fig. 6 Order assignment and phase interpolation for computer-simulated dense and sparse fringes. (a-1) and (a-2) Initial images; (b-1) and (b-2) Order assignment for skeletons; (c-1) and (c-2) Three-dimensional phase graphs by our method.

Download Full Size | PDF

The proposed methods not only obtain high accuracy of computer-simulated fringe pattern, but also give good results for experimentally-obtained optical fringe patterns with poor quality. Figure 7 shows two experimentally-obtained ESPI fringe pattern with poor quality. Figures 7(a-1) and 7(a-2) are the unloaded speckle pattern; Figs. 7(b-1) and 7(b-2) give the loaded speckle pattern; Figs. 7(c-1) and 7(c-2) show the ESPI fringe pattern with poor quality (which come from the subtraction of unloaded patterns and loaded patterns); Figs. 7(d-1) and 7(d-2) are the skeleton images corresponding to the fringe image of Figs. 7(c-1) and 7(c-2) where $Δ t = 0.2$ , $C = 0.0001$ for both the two groups of images, iterations are $n = 400$ and $n = 250$ , and the window size are 43 × 43 and 37 × 37, respectively; Figs. 7(e-1) and 7(e-2) show the results of phase interpolation by our method with $Δ t = 0.25$ for both the two groups of images, and $n = 6000$ , $n = 5000$ , respectively; Figs. 7(e-1) and 7(e-2) give the three-dimensional phase graph. To avoid the inherent edge distortion error and obtain the real skeleton pixels, we delete several white pixels at the boundary in Fig. 7(d-2). From Figs. 7(d-1) and 7(d-2), it can be seen that our proposed skeleton extraction method can work well on the optical fringe patterns with poor quality. Even if the original fringe image contrast is not obvious, our approach can still get the relative complete skeletons. The interpolation results shown in Figs. 7(e-1), 7(e-2), 7(f-1) and 7(f-2) illustrate that the proposed method can give the good phase information.

Fig. 7 skeleton extraction and phase interpolation for experimentally-obtained ESPI fringe pattern with poor quality. (a-1) and (a-2) The unloaded speckle pattern; (b-1) and (b-2) The loaded speckle pattern; (c-1) and (c-2) Fringe pattern with poor quality; (d-1) and (d-2) The skeleton images; (e-1) and (e-2) Gray-phase images by our method; (f-1) and (f-2) Three-dimensional phase graph.

Download Full Size | PDF

6. Conclusion and outlook

In this paper, we propose a skeleton extraction method and a phase interpolation method based on the use of PDEs. For skeleton extraction, the proposed method utilizes an anisotropic diffusion PDE to adjust the gradient vector field of the ESPI fringe image and judges the skeleton lines based on the divergence property of the adjusted gradient vector field. In addition, the proposed method uses the principle of heat conduction based on PDE to interpolate the whole-field phase. The main advantage of our method is that it can work on gray-scale ESPI fringe image directly without any preprocessing, thus avoiding information lost due to the inappropriate pre-processing, which is helpful to keep the continuity of the skeleton lines. The second advantage is the improved accuracy in phase extraction. The proposed method is able to yield whole-field phase for both the dense fringe pattern and the sparse fringe pattern. The third advantage is that the proposed phase interpolation method is faster than the method based on BP Neural Network and smoother than the RBF interpolation method.

Acknowledgment

This work is sponsored by National Nature Science Foundation of China (NSFC) under grant No. 61302127, Tianjin Science and Technology Supporting Projection under grant No. 13ZCZDGX02100 and No. 14ZCZDGX00033, Tianjin Research Program of Application Foundation and Advanced Technology under grant No.15JCYBJC16600, and Open Foundation of Key laboratory of Opto-electronic Information Technology of Ministry of Education (Tianjin University).

References and links

1. V. Bavigadda, R. Jallapuram, E. Mihaylova, and V. Toal, “Electronic speckle-pattern interferometer using holographic optical elements for vibration measurements,” Opt. Lett. 35(19), 3273–3275 (2010). [CrossRef] [PubMed]

2. G. Pedrini, W. Osten, and M. E. Gusev, “High-speed digital holographic interferometry for vibration measurement,” Appl. Opt. 45(15), 3456–3462 (2006). [CrossRef] [PubMed]

3. C. Quan, C. J. Tay, F. Yang, and X. He, “Phase extraction from a single fringe pattern based on guidance of an extreme map,” Appl. Opt. 44(23), 4814–4821 (2005). [CrossRef] [PubMed]

4. L. Lam, S. W. Lee, and C. Y. Suen, “Thinning methodologies-a comprehensive survey,” IEEE Trans. Pattern Anal. Mach. Intell. 14(9), 869–885 (1992). [CrossRef]

5. C. Direkoglu, Feature extraction via heat flow analogy (Dissertation, University of Southampton, 2009).

6. C. Tang, W. Lu, Y. Cai, L. Han, and G. Wang, “Nearly preprocessing-free method for skeletonization of gray-scale electronic speckle pattern interferometry fringe patterns via partial differential equations,” Opt. Lett. 33(2), 183–185 (2008). [CrossRef] [PubMed]

7. T. M. Lehmann, C. Gönner, and K. Spitzer, “Survey: interpolation methods in medical image processing,” IEEE Trans. Med. Imaging 18(11), 1049–1075 (1999). [CrossRef] [PubMed]

8. C. Tang, W. Lu, S. Chen, Z. Zhang, B. Li, W. Wang, and L. Han, “Denoising by coupled partial differential equations and extracting phase by backpropagation neural networks for electronic speckle pattern interferometry,” Appl. Opt. 46(30), 7475–7484 (2007). [CrossRef] [PubMed]

9. G. Wang, Y. J. Li, and H. C. Zhou, “Application of the radial basis function interpolation to phase extraction from a single electronic speckle pattern interferometric fringe,” Appl. Opt. 50(19), 3110–3117 (2011). [CrossRef] [PubMed]

10. Y. Wang, X. Ji, and Q. Dai, “Fourth-order oriented partial-differential equations for noise removal of two-photon fluorescence images,” Opt. Lett. 35(17), 2943–2945 (2010). [CrossRef] [PubMed]

11. F. Zhang, Z. Xiao, J. Wu, L. Geng, H. Li, J. Xi, and J. Wang, “Anisotropic coupled diffusion filter and binarization for the electronic speckle pattern interferometry fringes,” Opt. Express 20(20), 21905–21916 (2012). [CrossRef] [PubMed]

12. J. H. Jang and K. S. Hong, “A pseudo-distance map for the segmentation-free skeletonization of gray-scale images,” IEEE International Conference on Computer Vision, 2, 18–23 (2001).

13. Z. Yu and C. Bajaj, “A segmentation-free approach for skeletonization of gray-scale images via anisotropic vector diffusion,” IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 1, 415–420 (2004).

14. J. Zhou and J. Gu, “A model-based method for the computation of fingerprints’ orientation field,” IEEE Trans. Image Process. 13(6), 821–835 (2004). [CrossRef] [PubMed]

15. C. Tang, Z. Wang, L. Wang, J. Wu, T. Gao, and S. Yan, “Estimation of fringe orientation for optical fringe patterns with poor quality based on Fourier transform,” Appl. Opt. 49(4), 554–561 (2010). [CrossRef] [PubMed]

16. V. Espinosa-Duro, “Mathematical morphology approaches for fingerprint thinning,” International Carnahan Conference on Security Technology, 1, 43–45 (2002).

Skeleton extraction and phase interpolation for single ESPI fringe pattern based on the partial differential equations

Abstract

1. Introduction

2. Skeleton extraction for fringe patterns based on PDE

2.1 The initial gradient vector field

2.2 The diffusion of the gradient vector field

2.3 Skeleton extraction based on the divergence analysis

3. Phase interpolation for fringe patterns based on PDE

3.1 Assignment of the order number of skeletons

3.2 Phase interpolation

4. Numerical solution of PDEs

5. Experiments and results

6. Conclusion and outlook

Acknowledgment

References and links

Cited By

Figures (7)

Equations (21)

Optics Express