Separation of Two-Dimensional Mixed Circular Fringe Patterns Based on Spectral Projection Property in Fractional Fourier Transform Domain

: In this paper, a method for automatically separating the mixed circular fringe patterns based on the fractional Fourier transform (FrFT) analysis is proposed. Considering the mixed two-dimensional (2-D) Gaussian-based circular fringe patterns, detected by using an image sensor, we propose a method that can efﬁciently determine the number and parameters of each separated fringe patterns by using the FrFT due to the observed higher sparsity in the frequency domain than that in the spatial domain. First, we review the theory of FrFT and the properties of the 2-D circular fringe patterns. By searching the spectral intensities of the various fractional orders in the FrFT projected along both the frequency axes, the proposed method can automatically determine the total fringe number, the central position, binary phase, and the maximum fringe width of each 2-D circular fringe pattern. In the experimental results, both the computer-simulated and optically mixed fringe patterns are used to verify the proposed method. In addition, the additive Gaussian noise effects on the proposed method are investigated. The proposed method can still successfully separate the mixed fringe pattern when the signal-to-noise ratio is higher than 7 dB.


Introduction
In communication systems, when signals are disturbed by noise, we can use the Fourier transform (FT) method to transfer the signals to the frequency domain for noise filtering. However, in the case of mixed signals, their Fourier spectra may overlap, and thus we cannot directly separate the noise from the signals by using direct filtering in the frequency domain. In addition to the FT, the methods based on the fractional Fourier transform (FrFT) [1], in which the fractional orders can be sought for the best separation between the noise and signal spectra, have been found to be more advantageous than the FT in solving the spectral overlapping problems [2][3][4][5][6]. Previous studies have mostly focused on the separation of one-dimensional (time domain) mixed signals. Kutay and Ozaktas dealt with the problem of image restoration by using the two-dimensional (2-D) FrFT [7]. The FrFT of a 2-D signal f (x,y) can be defined as [8].
F α x α y [ f (x, y)](u, v) = K α x α y (x, y; u, v) f (x, y)dxdy and φ x = α x π/2 denotes the angle corresponding to the fractional order α x in the x direction. Note that the definition of K α y (y, u) is similar to that of K α x (x, u) and the range of α x and α y values is within [0,2].
Recently, optical technologies based on measuring the interference fringe patterns have been dramatically developed, especially for optical metrology, because the results are usually encoded as fringe patterns. In some cases, a single fringe pattern may contain multiple fringe sets. The fringe patterns are 2-D signals and are usually detected as 2-D images by using the image sensors. Chang et al. firstly use the FrFT to automatically separate the 2-D mixed Gaussian chirp signals [9]. Qian used the 2-D windowed Fourier transform to determine the phase and phase derivatives in fringe pattern analysis [10]. Trusiak et al. proposed several algorithmic solutions based on the notion of Hilbert-Huang transform for fringe pattern processing and analysis in optical interferometry [11]. Patorski [12] and Lai [13] proposed the methods for the mixed 2-D signals composed of the overlapping fringe patterns. Patorski's method utilizes the 2-D continuous wavelet transform (CWT) and is especially useful for separating the mixed Moiré patterns. However, improvement is still required for the cases of symmetrical circular fringe patterns. In Lai's method [13], the original mixed fringe patterns are firstly transformed to the spatial-frequency domain by using the 2-D FrFT. Then, the maxima projection method is proposed to search for the fractional Fourier spectra under all the possible fractional orders. The circular fringes are of the same width and their center positions can be determined according to the peak values detected in the spectra. With the determined fractional order, fringe width, and the center position, each of the separated circular fringes can be reconstructed.
This study deals with a more general case than the previous method for automatically separating the 2-D mixed multiple fringe patterns by using the FrFT. The previous method proposed by Lai [13] does not consider the condition of phase reversal phenomena among the fringe patterns. Thus, they only consider the fringe patterns, which are all in phase, and thus cannot reconstruct the mixed multiple fringe patterns, in which some are out of phase. In the proposed method, we extend the previous method to be able to deal with more general cases. A generalized formulation for 2-D circular fringe patterns is presented. However, only the fringe patterns with 0 and 180 degrees of phase differences are considered because, in some experimental cases, the fringe carrier is subject to contrast reversals. We determine the correlation coefficients between the separated signals with both cases and the original one. Therefore, the mixed 2-D fringe patterns with 0 and 180 degrees of phase differences can be successfully separated and reconstructed with acceptable quality.

Circular Fringe Pattern and Their Applications
The general formula of a 2-D fringe pattern with a constant phase shift ϕ 0 is shown as follows where I 0 denotes the background intensity, I 1 is the modulation amplitude, g(x, y) is a 2-D Gaussian distributed function and ϕ(x, y) denotes the phase function that can be used to generate a fringe pattern whose center is located at the coordinates (x 0 , y 0 ). Equations (4) and (5) show the definitions of ϕ(x, y) and g(x, y), respectively.
Here, a is the maximum width of the fringe pattern and σ is the standard deviation of the 2-D Gaussian function. Figure 1a,b show the two 2-D circular fringe patterns with 180 • phase difference in between, i.e., ϕ 0 = π. where I0 denotes the background intensity, I1 is the modulation amplitude, g(x, y) is a 2-D Gaussian distributed function and ( , ) denotes the phase function that can be used to generate a fringe pattern whose center is located at the coordinates (x0, y0). Equations (4) and (5) show the definitions of ( , ) and g(x, y), respectively.
Here, a is the maximum width of the fringe pattern and is the standard deviation of the 2-D Gaussian function. Figure 1a,b show the two 2-D circular fringe patterns with 180° phase difference in between, i.e., 0 = . The optimal fractional order αopt for a given Gaussian fringe pattern of size N by N can be determined as [10] (3), showing that the center of the fringe pattern is bright. On the other hand, Figure 1b demonstrates the case of 0 = in Equation (3), the center of the fringe pattern is dark and located at a different position. Figure 1c shows the mixed fringe pattern, which is the mixture of the ones shown in Figure  1a,b. Note that the contrast of the mixed fringe pattern has been decreased. If more circular fringe patterns are mixed, the number of fringes will not be easily recognizable.
There are various related works applying circular fringe patterns, especially in applications of the optical metrology. For low-quality fringe pattern enhancement and normalization, Trusiak et al. proposed an automatic selective reconstruction method to decompose a fringe pattern into a set of empirical modes using enhanced fast empirical mode decomposition [14]. Lu et al. proposed a method to estimate the parameters of optical fringes with a quadratic phase using the FrFT [15]. Sciammarella and Lamberti evaluated the basic assumptions of the standard methods utilized in fringe pattern analysis and presented the mathematical models to determine the displacement information encoded in fringe patterns [16]. Wu et al. proposed a method that combines the advantages of the FrFT with the least-squared fitting method in analyzing the fringe patterns [17]. Better parameter estimation can be achieved compared with the conventional FrFT method. To retrieve the phase information from a single-shot spatial carrier fringe pattern, Dong and Chen proposed an advanced FT analysis method [18], which can mitigate the problems in The optimal fractional order α opt for a given Gaussian fringe pattern of size N by N can be determined as [10] for continuous-time cases and for discrete-time cases, where f s is the sampling frequency in the 2-D spatial domain. Figure 1a demonstrates the case of ϕ 0 = 0 in Equation (3), showing that the center of the fringe pattern is bright. On the other hand, Figure 1b demonstrates the case of ϕ 0 = π in Equation (3), the center of the fringe pattern is dark and located at a different position. Figure 1c shows the mixed fringe pattern, which is the mixture of the ones shown in Figure 1a,b. Note that the contrast of the mixed fringe pattern has been decreased. If more circular fringe patterns are mixed, the number of fringes will not be easily recognizable.
There are various related works applying circular fringe patterns, especially in applications of the optical metrology. For low-quality fringe pattern enhancement and normalization, Trusiak et al. proposed an automatic selective reconstruction method to decompose a fringe pattern into a set of empirical modes using enhanced fast empirical mode decomposition [14]. Lu et al. proposed a method to estimate the parameters of optical fringes with a quadratic phase using the FrFT [15]. Sciammarella and Lamberti evaluated the basic assumptions of the standard methods utilized in fringe pattern analysis and presented the mathematical models to determine the displacement information encoded in fringe patterns [16]. Wu et al. proposed a method that combines the advantages of the FrFT with the least-squared fitting method in analyzing the fringe patterns [17]. Better parameter estimation can be achieved compared with the conventional FrFT method. To retrieve the phase information from a single-shot spatial carrier fringe pattern, Dong and Chen proposed an advanced FT analysis method [18], which can mitigate the problems in the traditional FT analysis method. To measure the whole-field out-of-plane deformation of targets, Ratnam et al. proposed a new class of circular fringe projection methods that employ circular fringes for three-dimensional shape measurement [19]. Consider the optical elements with full-field or circular pupils, and that the phase presents a quadratic behavior such as astigmatic, elliptical, or circular fringes. Muñoz-Maciel presented an FT-based method for phase recovery from a single interferogram with closed fringes [20]. Sciammarella et al. proposed a new method to extract the displacement information from fringe pattern in 2-D fields based on the property of FT or Hilbert transforms [21]. A method using the Simulated Annealing technique to determine the phase term from a single fringe pattern was proposed by Moré et al. [22]. In this method, the authors automatically partition the interferogram using a recursive method that stores in a quad tree data structure with a limit of fringes in each partition. Finally, Guo and Yang proposed a chirp FT-based method to estimate the parameters of Newton's rings [23]. In addition, the proposed method can also be used to measure the curvature radius of a spherical surface.
Recently, the applications of the methods for separating complex fringe patterns have been received a lot of attention. In addition to demodulating the phase of multiple fringe sets superimposed in a single image [11], for example, Pokorski and Patorski proposed the method based on the 2-D CWT to filter out parasitic fringes or to extract two separate superimposed fringe families with information encoded in their phase [24]. Trusiak and Patorski proposed an enhanced method for two-shot fringe pattern phase-amplitude demodulation using Gram-Schmidt orthonormalization with Hilbert-Huang transform by filtering out the spurious noise and background illumination and performing fringe normalization [25]. Another application is the method for nondestructive material evaluation by using two moiré interferometry fringe patterns with encoded orthogonal in-plane displacement information [26]. Therefore, aiming towards similar requirements, as shown in the above studies, we propose a method for automatically separating complex fringe patterns.

Proposed Method
According to the definition shown in Equation (1), the fractional orders α x and α y in the x and y directions, respectively, can be different. For simplicity, we assume that the fractional orders of Gaussian fringe patterns are identical in the proposed method. That is, Figure 2 shows the systematic block diagram of the proposed method.
The detailed procedures are listed as follows: 1. First, given an image of size N by N, which consists of mixed multiple 2-D fringe patterns, we apply the 2-D FrFTs to the whole image with the fractional orders of the range within [0.1] with the increment step 0.01 in either the x or y direction. When we perform the FrFTs in the x-direction, the order in the y-direction is fixed as 1; 2.
The orders with the maximum projection in the x-direction of the signal are calculated and recorded; 3.
Then, we perform the FrFTs in the y-direction and, similarly, the order in the xdirection is fixed as 1; 4.
The orders with the maximum projection in the y-direction of the signal are calculated and recorded as well; 5.
We calculate the orders of maximum projection, α x and α y , from the local maximum values of fractional Fourier spectra through a differential operation. Then, the set of best orders, α best , and the number of the paired orders in the set, M, can be determined by finding the intersections of the maximum projection orders α x,p and α y,p . Equations (8) and (9) show the definitions of the maximal projection orders for the x and y directions, respectively.
Here, we have assumed that the fractional orders of α x and α y should be identical in the FrFT. Therefore, the maximum width of for each fringe pattern, a, can be determined based on the best order α best by using Equation (6).

6.
By enumerating their absolute values and normalizing the transformed results F r (u,v), we deal with it by using a filter mask and searching the peak-value position in the filtered results. The filter mask is of size 3 × 3 and its weighting coefficients are The filtered image result is denoted as G(u,v) and its maximum g max is then determined. Through setting up a threshold value T between 0 and 1, we can get a binary map g bw from the filtered transformed result by using the criterion shown in Equation (10) In this binary map, we consider the neighboring area near the pixel value 1 as the same fringe pattern by searching the 8-connected positions. 7.
The original signal is FrFTed with the best order, α best , to calculate the center position of a specific circular fringe pattern. If the positions (u m , v m ) of the k f maximal values in the binary map derived from the filtered fractional Fourier spectra G(u,v) are detected, it means that we can separate k f fringe patterns from the original mixture signal; 8.
The center points (x 0 , y 0 ) in each corresponding circular fringe pattern can be determined according to Equation (11) 9.
After calculating the center points and maximum widths of the k f fringe patterns, we can reconstruct k f fringe patterns separately. For the cases of the fringe patterns with only 180 • phase difference, the reconstructed signals, s k1 , and s k2 (from the signal s k1 ) are generated by using Equation (3) with ϕ 0 = 0 o and ϕ 0 = π, respectively; 10. The correlation coefficients (CCs) between the original and reconstructed signals, s k1 , and s k2 , are denoted as cor k1 and cor k2 for all the k f reconstructed fringe patterns, respectively. The CC value r between the two 2-D fringe patterns X(i, j) and Y(i, j) of size N by N is defined in Equation (12), which is used to measure the similarity between the reconstructed and original signals.
Note that X and Y represent the average values of X(i, j) and Y(i, j), respectively. 11. By comparing the CC values of cor k1 and cor k2 , we can determine the constant phase shift ϕ 0 denoted in Equation (3). That is, we select the reconstructed signal that has a larger CC value than the other one; 12. Finally, all the k f separated fringe patterns can be solely reconstructed and the accuracy of each reconstructed fringe pattern can be evaluated by using the CC values.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 5 of 15 Step (6) FrFT Calculate CC value between s k1 and input signal, and so as s k2 , (cor k1 , cor k2 ) Step (5) Step (7) Calculate the center position (u m , v m ) of each circular fringe Step (8) Step (9) Step (10) Step (11) Step (1) Step (2) Step (3) Step (4) The detailed procedures are listed as follows: 1. First, given an image of size N by N, which consists of mixed multiple 2-D fringe patterns, we apply the 2-D FrFTs to the whole image with the fractional orders of the range within [0.1] with the increment step 0.01 in either the x or y direction. When we perform the FrFTs in the x-direction, the order in the y-direction is fixed as 1; 2. The orders with the maximum projection in the x-direction of the signal are calculated and recorded; 3. Then, we perform the FrFTs in the y-direction and, similarly, the order in the x-direction is fixed as 1; 4. The orders with the maximum projection in the y-direction of the signal are calculated and recorded as well; 5. We calculate the orders of maximum projection, αx and αy, from the local maximum values of fractional Fourier spectra through a differential operation. Then, the set of best orders, αbest, and the number of the paired orders in the set, M, can be determined by finding the intersections of the maximum projection orders ,p and ,p . Equations (8) and (9) show the definitions of the maximal projection orders for the x and y directions, respectively.

Experimental Results
The experiments are performed with MATLAB R2013b and a PC with Intel Core 2 Quad Q8200 CPU, DDR3 2G RAM, and Windows 7 OS. For simplicity, both the background intensity I 0 and the modulation amplitude I 1 in Equation (3) are set as unity. As shown in Figure 2, only two cases of the constant phase shifts, ϕ 0 = 0 • and ϕ 0 = 180 • are considered. Figure 3 shows the two mixed fringe patterns used in our experiments, which comprises six and eight 2-D fringe patterns, respectively. The image is of size 512 × 512. According to the block diagram shown in Figure 2, we apply the FrFTs on these pattern images with the fractional orders α = [0.01, 1] with an increment step 0.01 at the xand y-axis directions. The orders with the maximum projections in the xand y-direction of the signal are recorded as α x,p and α y,p , respectively.
Consider the mixed fringe pattern shown in Figure 3a. are considered. Figure 3 shows the two mixed fringe patterns used in our experiments, which comprises six and eight 2-D fringe patterns, respectively. The image is of size 512 × 512. According to the block d (a) (b) Consider the mixed fringe pattern shown in Figure 3a.    Consider the mixed fringe pattern shown in Figure 3a. Figure 4a-f shows part of the projection results with the various fractional orders. Figure 4a-c shows the projection results in the x-direction when the order is 1 in the y-direction, while Figure 4d-f shows the results in the y-direction when the order is 1 in the x-direction. By applying the derivative operations on the projection values shown in Figure 5, the fractional orders of the local maxima in the x and y directions ar  To find out how many possible fringe patterns are contained in the mixed one, we determine the fractional Fourier spectra (i.e., intensities of the FrFT) of the signals with the determined best orders and corresponding fringe widths. By applying a peak detection filter, we find six peaks gmax and the corresponding positions of fringe patterns, gbw in the four spectra. The six separated fringe patters can thus be reconstructed. Figure 6a,c,e,g show the normalized fractional Fourier spectra G(u,v) for the determined best fractional orders αbest = {0.73, 0.76, 0.79, 0.82} and the filtered results are shown in Figure 6b,d,f,h. determined best orders and corresponding fringe widths. By applying a peak detection filter, we find six peaks g max and the corresponding positions of fringe patterns, g bw in the four spectra. The six separated fringe patters can thus be reconstructed. Figure 6a,c,e,g show the normalized fractional Fourier spectra G(u,v) for the determined best fractional orders α best = {0.73, 0.76, 0.79, 0.82} and the filtered results are shown in Figure 6b,d,f,h. There are six peak values observed in the filter results, which correspond to six single fringe patterns in the mixed one. By using Equations (9) and (10), we can determine the maxima in the filter results and then the center positions. With the maximum widths of fringe patterns, a, the six separated circular fringe patterns can be reconstructed and shown in Figure 7a-f. Figure 8a-h shows the eight separated circular fringe patterns from a more complicated mixed signal shown in Figure 3b, based on the similar process shown in Figures 4-6.
Finally, the results are shown in Tables 1-3. Figure 7a-f shows the separated fringe patterns, which are reconstructed by using the determined parameters shown in Tables 1-3. In Table 3, we can find that the CC values for most reconstructed signals are higher than 0.746. However, the CC value r of the reconstructed signal in Figure 7f is only 0.431, which is much lower than the others and could be caused by the higher variability on the determined position and width. Tables 4 and 5 show the center positions and the maximal widths, respectively, of the eight circular fringe patterns separated from Figure 3b. As shown in Table 4, the differences between the original and determined center positions are small and the range is similar to that shown in Table 2. Regarding the difference of the maximal width of the circular fringe patterns shown in Table 5, the error percentages are all less than 3%. The CC values between the eight reconstructed and the original circular fringe patterns are also shown in Table 5, which are, on average, lower than that shown in Table 3. The possible reasons for this are the increased number of mixed circular fringe patterns and the partial overlapping among these patterns. Table 1. The determined best orders α best and width a of fringe patterns.

Parameter
Result      Tables 1-3. Figure 7a-f shows the separated fringe patterns, which are reconstructed by using the determined parameters shown in Tables 1-3. In Table 3, we can find that the CC values for most reconstructed signals are higher than 0.746. However, the CC value r of the reconstructed signal in Figure 7f is only 0.431, which is much lower than the others and could be caused by the higher variability on the determined position and width. Tables 4 and 5 show the center positions and the maximal widths, respectively, of the eight circular fringe patterns separated from Figure 3b. As shown in Table 4, the differences between the original and determined center positions are small and the range is similar to that shown in Table 2. Regarding the difference of the maximal width of the circular fringe patterns shown in Table 5, the error percentages are all less than 3%. The CC values between the eight reconstructed and the original circular fringe patterns are also shown in Table 5, which are, on average, lower than that shown in Table 3. The possible reasons for this are the increased number of mixed circular fringe patterns and the partial overlapping among these patterns.  Figure 9a shows a multiple-fringe pattern, which is the recorded intensity of light field diffracted from three pinholes illuminated by a plane wave in an optical experiment. The wavelength of blue laser illumination is 473 nm. Figure 9b shows the grayscale image of Figure 9a. By using the proposed method, three fringe patterns can be separated from the mixed signal. Figure 10a-c shows the separated circular fringes by using the proposed method, while Figure 10d shows the mixed fringe pattern reconstructed by combining the three separated ones. This example shows that the proposed method can be applied to the separation problem for the real fringe patterns. The noise effects on the mixed fringe patterns in the proposed method are also investigated. Here we consider the additive Gaussian noise with various power. Let the signal to noise ratio (SNR) be defined as SNR (in dB) = 10 log 10 P s P n (13) where P s and P n denote the signal power and noise power, respectively. Figure 11a-d show the noisy fringe patterns, which are the original fringe pattern in Figure 3a with additive Gaussian noise of SNR values 1, 3, 5, and 7 dB, respectively. Then, we test the proposed method with noisy fringe patterns of the SNR values from 1 to 10 dB. Figure 9a shows a multiple-fringe pattern, which is the recorded intensity of light field diffracted from three pinholes illuminated by a plane wave in an optical experiment. The wavelength of bl The noise effects on the mixed fringe patterns in the proposed method are also investigated. Here we consider the additive Gaussian noise with various power. Let the signal to noise ratio (SNR) be defined as where Ps and Pn denote the signal power and noise power, respectively. Figure 11a-d show the noisy fringe patterns, which are the original fringe pattern in Figure 3a with additive Gaussian noise of SNR values 1, 3, 5, and 7 dB, respectively. Then, we test the proposed method with noisy fringe patterns of the SNR values from 1 to 10 dB.  Figure 9a shows a multiple-fringe pattern, which is the recorded intensity of light field diffracted from three pinholes illuminated by a plane wave in an optical experiment. The wavelength of blue laser illumination is 473 nm. Figure 9b shows the grayscale image of Figure 9a. By using the proposed method, three fringe patterns can be separated from the mixed signal. Figure 10a-c shows the separated circular fringes by using the proposed method, while Figure 10d shows the mixed fringe pattern reconstructed by combining the three separated ones. This example shows that the proposed method can be applied to the separation problem for the real fringe patterns. The noise effects on the mixed fringe patterns in the proposed method are also investigated. Here we consider the additive Gaussian noise with various power. Let the signal to noise ratio (SNR) be defined as where Ps and Pn denote the signal power and noise power, respectively. Figure 11a-d show the noisy fringe patterns, which are the original fringe pattern in Figure 3a with additive Gaussian noise of SNR values 1, 3, 5, and 7 dB, respectively. Then, we test the proposed method with noisy fringe patterns of the SNR values from 1 to 10 dB. The noise effects on the mixed fringe patterns in the proposed method are also investigated. Here we consider the additive Gaussian noise with various power. Let the signal to noise ratio (SNR) be defined as where Ps and Pn denote the signal power and noise power, respectively. Figure 11a-d show the noisy fringe patterns, which are the original fringe pattern in Figure 3a with additive Gaussian noise of SNR values 1, 3, 5, and 7 dB, respectively. Then, we test the proposed method with noisy fringe patterns of the SNR values from 1 to 10 dB.  Table 6 shows the detected numbers of separated fringes and the corresponding CC values between the original and reconstructed fringe patterns under different SNR values of the noisy mixed fringe patterns. As shown in Table 6, the number of detected fringe patterns varies when the SNR values are within 1 and 7 dB. When the SNR values are 1, 2, 3, and 6 dB, some of the fringe patterns are not detected and are denoted as not available (NA) in the tables. On the other hand, some extra fringe pattern can be detected when the SNR value is 4 dB. When the SNR value is 5 dB, the CC values the fringe patterns No. 4 and No. 5 are −0.009 and 0.004, respectively, which means that these two fringe patterns are not missing. Instead, two extra fringe patterns, at different locations, are detected. When the SNR value is 7 dB, the fringe number and the CC values are identical to the noiseless case. Note that the test results (both the fringe number and CC values) for SNR values higher than 7 dB are the same.

Conclusions
We propose an automatic signal separation method for separating the 2-D mixed circular fringe patterns and considering the possible phase reversal condition. The proposed method utilizes the fractional Fourier transform to search for the best fractional orders so that the peak values corresponding to the single circular fringe pattern can be observed in the spectral domain. A filtering scheme is used to determine the precise center position of each fringe pattern. Then, the separated fringe pattern can be reconstructed based on the determined center position and maximum width. We have also demonstrated that both the digital and optical mixed fringe patterns can be successfully separated and most of the reconstructed signals have high correlation coefficients compared to the original ones. In addition, the proposed method can successfully separate the mixed fringe patterns with additive Gaussian noise when the SNR values is as high as 7 dB.
In our future work, the sampling rate effects, and the minimum detectable fringe size of the proposed method, will be investigated. In addition, we will extend the capability of proposed method, so that the fringe pattern with arbitrary constant phase shift ϕ 0 other than 180 • can be detected as well.