Defect detection using windowed Fourier spectrum analysis in diffraction phase microscopy

Automated defect identification is an important problem in non-destructive testing and evaluation, and is commonly needed for applications such as experimental mechanics, material inspection and condition-based maintenance. Optical interferometric techniques such as electronic speckle pattern interferometry [1], digital holographic interferometry [2] and shearography [3] are popular for non-contact and high precision detection of defects. In these techniques, we obtain a fringe pattern or interferogram which encodes information about the quantity of interest. The main objective for defect detection algorithms is to identify defects such as the presence of cracks, fracture and debonds from a fringe pattern [4, 5]. The fringe pattern serves as a good marker for defect identification because of the rapid fringe density variations in the vicinity of defect region. For identifying defects from fringe patterns, fringe tracking and contouring based methods [6, 7] were proposed, but they involve tedious and error-prone scanning and iterative operations. Another popular approach has been to use pattern classification [8–10] and neural networks for fault detection [11]. However, these methods require extensive training sets and prior data obtained from fringe patterns with different values of defect, noise and loading parameters which limits the performance and practical feasibility of these methods. Phase-shifting based approaches [12, 13] for defect identification have been reported but they require several phase-shifted fringe patterns to be captured which can be tedious experimentally. Similarly, methods based on numerical phase differentiation [14] and Fourier transform [15] have been proposed but these operations are usually sensitive to noise.

Another class of techniques for defect detection from fringe patterns include the spatial frequency analysis methods based on windowed Fourier transform (WFT) [16], wavelet transform (WT) [17] and Wigner-Ville distribution [18]. The windowed Fourier transform method relies on normalized cross-correlation approach where a database of WFT elements is constructed, and these elements are correlated with a reference and test fringe patterns. The wavelet transform method relies on computing the correlation between the fringe pattern and a mother wavelet at a given scale value. On the other hand, the Wigner-Ville distribution based method relies on computing a normalized measure of spatial phase derivatives from the reference and test fringe patterns. A major limitation of these methods is the requirement of a predetermined or preset threshold value to highlight the defect region which requires careful control and monitoring and is operator dependent. In addition, the procedure involved is usually computationally intensive.

In this paper, we propose an elegant method for defect detection which eliminates the requirement of multiple phase-shifted fringe patterns, multiple training data sets or prior knowledge of thresholds, and is also computationally efficient and robust against noise. The effectiveness and utility of the proposed technique is demonstrated in diffraction phase microscopy [19, 20] which is a common-path interferometric technique for non-invasive metrology. The theory of the proposed method is outlined in section 2. The simulation and experimental results are demonstrated in section 3 followed by discussions and conclusions.

Consider a carrier fringe pattern which can be mathematically expressed as [21]

$$\begin{eqnarray}&&I(x,y)=a(x,y)+b(x,y)\cos [\phi (x,y)+{\omega }_{{cx}}x+{\omega }_{{cy}}y]\end{eqnarray} \tag{ 1 }$$

where I(x, y) is the recorded intensity, a(x, y) is the background intensity, b(x, y) is the fringe amplitude and ϕ(x, y) is the phase distribution. Here, ω_cx and ω_cy are angular carrier frequencies along x and y directions. Carrier fringe patterns with mathematical form described above are frequently obtained in optical techniques such as fringe projection [22, 23], off-axis interferometry [2] and diffraction phase microscopy [20]. The first step in our method is to retrieve the complex analytic signal from the fringe pattern by computing a two-dimensional Fourier transform of the fringe pattern and applying spectral filtering to preserve only the side lobe centered around the carrier frequency [22]. The carrier frequency can be estimated by tracing the peak of the Fourier spectrum. After amplitude normalization, this yields a complex signal of the form

$$\begin{eqnarray}&&f(x,y)=\exp [j(\phi (x,y)+{\omega }_{{cx}}x+{\omega }_{{cy}}y)]\end{eqnarray} \tag{ 2 }$$

In the next step, we compute the windowed Fourier transform of the above signal as [24–26]

$$\begin{eqnarray}\begin{array}{rcl}{Wf}(u,v,\xi ,\eta ) & = & \iint f(x,y)g(x-u,y-v)\\ & & \exp [-j(\xi x+\eta y)]\,{dxdy}\end{array}\end{eqnarray} \tag{ 3 }$$

In the above equation, Wf(u, v, ξ, η) is a four dimensional signal having u and v as the spatial coordinates, and ξ and η as the angular frequency coordinates. Further, g(x, y) is the two-dimensional window function, which effectively enables local processing of the fringe pattern. Substituting the value of f(x, y) in above equation, we have

$$\begin{eqnarray}\begin{array}{rcl}{Wf}(u,v,\xi ,\eta ) & = & \iint \exp [j(\phi (x,y)+{\omega }_{{cx}}x+{\omega }_{{cy}}y)]\\ & & g(x-u,y-v)\exp [-j(\xi x+\eta y)]\,{dxdy}\end{array}\end{eqnarray} \tag{ 4 }$$

Using the transformations τ₁ = (x − u) and ${\tau }_{2}=(y-v)$, the above equation is modified as,

$$\begin{eqnarray}\begin{array}{rcl}{Wf}(u,v,\xi ,\eta ) & = & \iint \exp [j\phi (u+{\tau }_{1},v+{\tau }_{2})]\\ & & \exp [j({\omega }_{{cx}}(u+{\tau }_{1})+{\omega }_{{cy}}(v+{\tau }_{2}))]g({\tau }_{1},{\tau }_{2})\\ & & \exp [-j(\xi (u+{\tau }_{1})+\eta (v+{\tau }_{2}))]\,d{\tau }_{1}d{\tau }_{2}\end{array}\end{eqnarray} \tag{ 5 }$$

Using first order Taylor series expansion of phase within the window region, we have,

$$\begin{eqnarray}&&\phi (u+{\tau }_{1},v+{\tau }_{2})=\phi (u,v)+{\tau }_{1}{\phi }_{u}(u,v)+{\tau }_{2}{\phi }_{v}(u,v)\end{eqnarray} \tag{ 6 }$$

where ϕ_u(u, v) = ∂ ϕ(u, v)/∂ u and ϕ_v(u, v) = ∂ ϕ(u, v)/∂ v are the spatial phase derivatives. Substituting the above expression, we have the windowed Fourier transform as,

$$\begin{eqnarray}\begin{array}{rcl}{Wf}(u,v,\xi ,\eta ) & = & {e}^{j(\phi (u,v)+({\omega }_{{cx}}-\xi )u+({\omega }_{{cy}}-\eta )v)}\\ & & \iint g({\tau }_{1},{\tau }_{2}){e}^{j({\tau }_{1}{\phi }_{u}(u,v)+{\tau }_{2}{\phi }_{v}(u,v))}\\ & & {e}^{-j[(\xi -{\omega }_{{cx}}){\tau }_{1}+(\eta -{\omega }_{{cy}}){\tau }_{2}]}\,d{\tau }_{1}d{\tau }_{2}\end{array}\end{eqnarray} \tag{ 7 }$$

which further reduces to

$$\begin{eqnarray}\begin{array}{rcl}{Wf}(u,v,\xi ,\eta ) & = & {e}^{j(\phi (u,v)+({\omega }_{{cx}}-\xi )u+({\omega }_{{cy}}-\eta )v)}\\ & & \iint g({\tau }_{1},{\tau }_{2}){e}^{-j[(\xi -{\omega }_{{cx}}-{\phi }_{u}(u,v)){\tau }_{1}+(\eta -{\omega }_{{cy}}-{\phi }_{v}(u,v)){\tau }_{2}]}d{\tau }_{1}d{\tau }_{2}\end{array}\end{eqnarray} \tag{ 8 }$$

Denoting the two-dimensional (2D) Fourier transform of the window function as G(ω_x, ω_y), the term inside the double integral in the above equation is effectively the 2D Fourier transform of the window function g(τ₁, τ₂) evaluated at the angular spatial frequencies ${\omega }_{x}=\xi -{\omega }_{{cx}}-{\phi }_{u}(u,v)$ and ${\omega }_{y}=\eta -{\omega }_{{cy}}-{\phi }_{v}(u,v)$. Accordingly, the windowed Fourier spectrum can be expressed as

$$\begin{eqnarray}&&| {Wf}(u,v,\xi ,\eta )| =| G(\xi -{\omega }_{{cx}}-{\phi }_{u}(u,v),\eta -{\omega }_{{cy}}-{\phi }_{v}(u,v))| \end{eqnarray} \tag{ 9 }$$

For our analysis, a Gaussian window of the following form was chosen,

$$\begin{eqnarray}&&g(x,y)=\exp \left[-\left(\displaystyle \frac{{x}^{2}}{2{\sigma }_{x}^{2}}+\displaystyle \frac{{y}^{2}}{2{\sigma }_{y}^{2}}\right)\right]\end{eqnarray} \tag{ 10 }$$

where σ_x and σ_y are parameters related to the size of the window L_x (along x) and L_y (along y) such that ${L}_{x}=6{\sigma }_{x}+1$ and ${L}_{y}=6{\sigma }_{y}+1$.

Next, we evaluate the windowed Fourier spectrum at the carrier frequency by substituting ξ = ω_cx, η = ω_cy, which results in

$$\begin{eqnarray}&&| {Wf}(u,v,{\omega }_{{cx}},{\omega }_{{cy}})| =| G(-{\phi }_{u}(u,v),-{\phi }_{v}(u,v))| \end{eqnarray} \tag{ 11 }$$

The above equation is the key mathematical formulation for our defect detection method. As the window g(x, y) is of Gaussian form, its two-dimensional Fourier transform G(ω_x, ω_y) is also a Gaussian function. Importantly, the Gaussian window has a low pass behavior such that $| G(0,0)|$ has the peak spectral value, and value of the function $| G({\omega }_{x},{\omega }_{y})|$ decreases rapidly for spatial frequencies farther from the zero frequency. For a defect containing fringe pattern, the defect region is characterized by high fringe density or equivalently high spatial frequencies [18]. Hence, the phase derivatives ϕ_u(u, v) and ϕ_v(u, v) are markedly larger in the region associated with defects as compared to a defect-free region. Consequently, by considering the low pass behavior of the Gaussian window and the high spatial phase derivative value in the vicinity of defect, we have the following inequality

$$\begin{eqnarray}&&\mathop{\underbrace{| {Wf}(u,v,{\omega }_{{cx}},{\omega }_{{cy}})| }}\limits_{u,v\,\in \,{\bf{defect}}}\lt \mathop{\underbrace{| {Wf}(u,v,{\omega }_{{cx}},{\omega }_{{cy}})| }}\limits_{u,v\,\not\hspace{-2pt}{\in }{\bf{defect}}}\end{eqnarray} \tag{ 12 }$$

In other words, the windowed Fourier spectrum evaluated at the carrier frequency is significantly lower in the region with defect as opposed to a defect-free region. In the proposed method, this property is utilized for defect detection.

In the next step, we applied an automated thresholding criterion for defect detection from the aforementioned windowed Fourier spectrum based on Rosin's histogram algorithm [27]. For a fringe pattern with defects, we assume that the defect containing region is relatively smaller in area compared to the defect-free or background region, which is usually true for many practical cases. Consequently, the histogram of the computed WFT spectrum is predominantly unimodal, and an automated unimodal thresholding criterion such as the Rosin's algorithm is preferable.

To gain insights about the working of the proposed method, consider a simulated fringe pattern (size 1024 × 1024 pixels) with 'eye' type defect as shown in figure 1(a). The fringe pattern is corrupted with additive white Gaussian noise at signal to noise ratio (SNR) of 5 dB. The windowed Fourier spectrum was evaluated at the carrier frequency and normalized (peak value scaled to unity) for ease of visualization. The normalized spectrum is shown in figure 1(b). The histogram of the normalized spectrum is shown in figure 1(c). For better understanding, the semilog plot of the histogram is also shown in figure 1(d). The histogram serves as the input to Rosin's thresholding algorithm. In the algorithm, a line segment is drawn from the peak of the histogram curve to the first minimum. We denote the minimum point on the histogram curve as p₁ with coordinates (x₁, y₁), and the peak as point p₂ with coordinates (x₂, y₂). The equation of the line joining these two points can be modeled as $y={a}_{\mathrm{1,2}}x+{b}_{\mathrm{1,2}}$ where the slope ${a}_{\mathrm{1,2}}=({y}_{2}-{y}_{1})/({x}_{2}-{x}_{1})$ and the intercept ${b}_{\mathrm{1,2}}={y}_{1}-{a}_{\mathrm{1,2}}{x}_{1}$. Subsequently, a point in the histogram curve is located where the perpendicular distance from the line segment is maximum. This point constitutes the threshold in Rosin's algorithm. Mathematically, the length of a line segment drawn from an arbitrary point on the histogram with coordinates (x₀, y₀) to the straight line described above is given as

$$\begin{eqnarray}&&d({x}_{0},{y}_{0})=\displaystyle \frac{| {y}_{0}-{a}_{\mathrm{1,2}}{x}_{0}-{b}_{\mathrm{1,2}}| }{\sqrt{1+{a}_{1,2}^{2}}}\end{eqnarray} \tag{ 13 }$$

In Rosin's algorithm, the threshold is the point T at which d in the above equation is maximum with respect to the coordinates $({x}_{0},{y}_{0})$ on the curve. Numerically, the maximum length of the perpendicular line from the line segment is calculated using an iterative method [28].

Zoom In Zoom Out Reset image size

After computing the threshold value, the intensity of each pixel present in the normalized windowed Fourier spectrum image is compared with the threshold. In our method, the pixels with intensity values below the threshold are considered to be associated with defects, whereas the remaining pixels are assumed to be related to the background. Further, all pixels with spectrum values below the threshold are assigned the value of unity and rest are assigned the value of zero, effectively leading to a binary image. This binary conversion of the normalized windowed Fourier spectrum leads to a better visualization of the defect affected region which is clearly evident in figure 1(e).

In figure 2(a), we present a fringe pattern without defects simulated with SNR of 5 dB. The normalized windowed Fourier spectrum for this signal is shown in figure 2(b). From the figure, it is evident that in the absence of a defect, the windowed Fourier spectrum is essentially uniform with little variation. The corresponding histogram along with Rosin's threshold (marked T) is shown in figure 2(c) and the corresponding semilog plot is shown in figure 2(d). The resulting binary image is shown in figure 2(e).

Zoom In Zoom Out Reset image size

Note that due to the finite size of the fringe pattern and non-availability of data beyond the corners, the windowing operation leads to errors at the borders of the fringe pattern. Hence, for all subsequent figures related to our analysis, the border pixels are neglected to ignore the errors at the boundaries.

3.1. Simulation

The performance of the proposed method is tested for identifying four different types of faults or defects: 'eye', 'groove', 'bend' and 'compression', which occur frequently in fringe patterns [16]. Further, additive white Gaussian noise is added to the fringe pattern to validate the robustness of the method against noise. To obtain a quantifiable measure for the method's performance, we also performed statistical analysis by computing the true positive rate (TPR), also referred to as sensitivity, and true negative rate (TNR), also referred to as specificity, which are standard measures of quantifying success for defect identification. TPR signifies what proportion of defect is detected, and TNR indicates how specifically the defect is detected. Mathematically, TPR and TNR are expressed as

$$\begin{eqnarray}&&{TPR}=\displaystyle \frac{{TP}}{{TP}+{FN}}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{TNR}=\displaystyle \frac{{TN}}{{TN}+{FP}}\end{eqnarray} \tag{ 15 }$$

where TP indicates true positives or number of pixels correctly identified as pertaining to defects, TN indicates true negatives or number of pixels correctly identified as not pertaining to defects, FP indicates false positives or number of pixels wrongly classified as pertaining to defects, and FN indicates false negatives or number of pixels wrongly classified as not pertaining to defects, [29]. Both TPR and TNR values range between 0 and 1. In addition, both TPR and TNR should be close to unity for successful detection of defect.

For comparison, we also performed defect identification from the simulated fringe patterns using the standard wavelet transform method [17]. For this method, we evaluated the wavelet transform of the defect containing fringe pattern by using the Mexican hat mother wavelet at a scale value of unity and zero rotation angle. In the wavelet transform method, the spectrum values are higher in the vicinity of defect. Hence, in the normalized wavelet spectrum image, pixels having intensity higher than a prior determined threshold are assigned a value of unity, whereas the remaining pixels are allotted zero values which leads to a binary image. As the method requires a operator supplied threshold, we applied the wavelet transform method at three different threshold values denoted by $T1,T2$, and T3.

All simulated fringe patterns have size 1024 × 1024 pixels. For the proposed method, we used a Gaussian window of size 91 × 91 pixels. All computations were performed using Python programming language on a workstation with Intel i5 processor and 16 GB memory. In all cases, the binary image serves as the marker for fault identification since the pixels related to defect regions are shown has having high or unity value whereas the defect-free background is shown as having low or zero value. The binary image is also utilized for computing the TPR and TNR values.

The noisy fringe pattern simulated with SNR of 3 dB and 'eye' pattern fault is shown in figure 3(a). The normalized windowed Fourier spectrum computed by the proposed approach is shown in figure 3(b). The binary image obtained after applying Rosin's automated thresholding algorithm is shown in figure 3(c). This figure clearly highlights the defect detection capability of the proposed method even in the presence of high noise. For comparison, the binary images obtained from the wavelet transform method using threshold values of T1 = 50%, T2 = 60% and T3 = 70% of the peak wavelet spectrum value are shown in figures 3(d)–(f). Depending on the choice of the threshold value, the wavelet transform method exhibits varying degrees of defect identification performance, which is clear from these binary images. Further, the TPR and TNR curves plotted with respect to SNRs ranging from 0 dB to 20 dB for the proposed method, and the wavelet transform method computed with different threshold values is shown in figures 3(g) and (h). The proposed method shows a uniformly high value for both TPR and TNR even with varying noise, whereas the wavelet transform method exhibits a TPR and TNR trend dependent on both the noise level and selected threshold. Importantly, the wavelet transform method shows poor values of TPR and TNR especially for low SNRs.

Zoom In Zoom Out Reset image size

The performance of the proposed method and wavelet transform method for fringe patterns with 'groove' type defect are shown in figure 4. Similarly, the results for identifying 'bend' and 'compression' type defects from fringe patterns are shown in figures 5 and 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is clear from the above figures that the proposed method offers superior performance over the standard wavelet transform method for identifying defects from fringe patterns with respect to robustness against noise. This is further evident from the uniformly close to unity values of TPR and TNR curves obtained with the proposed method at different noise levels. In contrast, the performance of the wavelet transform method is significantly deteriorated with noise. Moreover, the effectiveness of the wavelet transform method is dependent on the selected threshold value, which is a major drawback of all methods requiring preset thresholds.

The main limitation of the proposed method is identifying small defects at severe levels of noise. To illustrate this, noisy fringe patterns with small and moderately large eye defects were simulated at SNRs of −5 dB, 0 dB and 5 dB, which are shown in figures 7(a)–(c). For better clarity, the regions containing the small defect in the fringe patterns and marked by red borders are cropped and highlighted in figures 7(d)–(f). The normalized windowed Fourier spectra for the three cases are shown in figures 7(g)–(i). Finally, the defects identified using the proposed method are shown in figures 7(j)–(l). For the cases of high noise levels, i.e. SNRs at −5 dB and 0 dB, it is evident from figures 7(j) and (k) that the small defect is not clearly identified by the proposed method. Further, some regions were wrongly identified as defect containing regions in these figures leading to false positives. On the other hand, at SNR of 5 dB, the small defect is identified by the proposed method. In addition, the moderately large defect is correctly identified in all three cases. This shows that severe noise can deteriorate the sensitivity of the proposed method especially for identifying small defects in fringe patterns.

Zoom In Zoom Out Reset image size

3.2. Experiment

We also tested the practical validity of the proposed method in diffraction phase microscopy, which relies on a common-path optical interferometric setup to minimize the detrimental effects of vibrations. The schematic of the system is shown in figure 8. In diffraction phase microscopy, coherent light is incident on the sample under test and the wave scattered from the sample is imaged via an objective lens (OL) and tube lens (TL) assembly on a diffraction grating (G). The grating creates multiple diffraction orders, among which only the zeroth and first order waves are allowed to pass further; other orders are blocked. Both these waves pass through a 4f lens system (L1,L2). The first order wave is filtered using a pinhole (P) in the Fourier plane of first lens, so that it resembles a plane wave after emerging from the second lens. In contrast, the zero order wave passes unfiltered or unmodified. The camera records the interference of these waves and provides the interferogram or fringe pattern. Detailed implementation of the experimental setup is discussed in [30]. We applied the proposed technique for identifying faults associated with sample preparation and processing processes such as photolithography and thermal evaporation. For each case, we recorded the defect containing interferogram and a reference interferogram (without defects) which is used to minimize the effects of system aberrations and unwanted diffraction patterns induced by dust particles etc. These experimental abnormalities are common to both the interferograms, and thus can be effectively eliminated.

Zoom In Zoom Out Reset image size

For sample preparation, a thin film of Aluminium of thickness 150 nm was initially coated over silicon wafer substrate. A photoresist was spin coated over Aluminium and a specific area of the photoresist was exposed using photolithography. Subsequently, the Aluminium and remaining photoresist were removed using etchant and developer solution. This process created a surface pattern or defect, which is captured by the interferogram recorded in diffraction phase microscopy, and is shown in figure 9(a). In figure 9(b), the corresponding reference interferogram obtained by imaging defect-free region of the sample is shown. The normalized windowed Fourier spectrum is shown in figure 9(c). Finally, the detected defect using the proposed method is highlighted in figure 9(d).

Zoom In Zoom Out Reset image size

We also considered another sample for testing where an Aluminium thin film of thickness 200 nm is coated over glass substrate using thermal evaporation. Thermal evaporation is a widely used technique for thin film coating. However, there may be areas where the coating material is not properly deposited on the substrate which leads to defects. Such defects are visible in the interferogram shown in figure 10(a). In this figure, the interferogram clearly shows an area with improper, non-uniform and incomplete coating, and thus emerges as a mostly dark region since the incident light is simply transmitted through the glass substrate with minimal reflection. The reference interferogram is shown in figure 10(b). The normalized windowed Fourier spectrum is shown in figure 10(c) and the detected defect by the proposed method is shown in figure 10(d).

Zoom In Zoom Out Reset image size

These results clearly validate the utility of the proposed method for automated and microscopic defect detection in diffraction phase microscopy.

The proposed method offers a fast and elegant approach for defect detection in optical metrology. Compared to the method proposed by Kemao [16], our method does not rely on computing correlation between multiple WFT elements between a reference and test fringe patterns or the WFT ridges which are computationally intensive operations. As a result, the proposed method is highly computationally efficient, and a 1024 × 1024 interferogram was processed within two seconds. Further, in contrast to the pattern classification and phase-shifting methods, the proposed method eliminates the need for capturing multiple fringe patterns which leads to high throughput. The proposed method also exhibits high robustness against noise, which was extensively tested using the TPR and TNR based analysis. Note that for a large window size, more data samples are captured for processing which leads to better performance against noise but the linear phase approximation would be unsuitable. On the other hand, an extremely small window might lead to high sensitivity towards noise. Hence, there is great scope for exploring adaptive windowing scheme in this domain. In addition, though we used spectral peak tracing method for carrier detection, other carrier retrieval techniques [31] can also be explored. Finally, the proposed method does not require preset thresholds for defect identification, unlike existing space-frequency methods, which leads to operational ease and simplicity for the user.

The authors proposed a fast, robust and high throughput method for automated defect detection, and demonstrated its application in diffraction phase microscopy. The method's utility was validated using simulation and experimental results. The authors believe that the method has immense potential for optical interferometry based non-destructive testing and inspection.

The work was supported by funding received from ISRO-IITK-STC under grant number STC0146.