Continuous phase estimation from noisy fringe patterns based on the implicit smoothing splines

Maciek Wielgus; Krzysztof Patorski; Pablo Etchepareborda; Alejandro Federico

doi:10.1364/OE.22.010775

1. Introduction

In optical metrological systems the measured quantity modifies such light beam parameters as its amplitude (intensity), propagation direction, frequency, phase and polarization state. Decoded light intensity changes are related to the measurand. In the case of studying phase modulating objects (transmitting or reflecting) interferometry constitutes well established and powerful indirect demodulation technique frequently used in science and technology. The so-called full-field measurement techniques with parallel data acquisition and processing for all objects points imaged onto the 2D matrix sensor deliver the output in the form of a fringe pattern, e.g., an interferogram (i.e., classical two-beam, holographic, speckle correlogram), moiregram, elastooptic polariscope image, and raster or structured illumination. Note that each fringe pattern has its characteristic features strongly dependent on the optical method used (e.g., coherent or noncoherent) and object under test (e.g., with specular of diffuse reflection).

The automated fringe pattern analysis (AFPA) enables optical methods to be user friendly, reliable and accurate [1, 2]. For out-of-laboratory tests single frame techniques are highly desired. They can be implemented using much simpler experimental hardware, are much more robust with respect to environmental disturbances and facilitate investigation of transient events. Simplified setup configurations call for more sophisticated software (algorithmic) solutions to provide demanded measurement accuracy. Among single frame AFPA techniques being dynamically developed over the last two decades we have: the Fourier (FT) and windowed Fourier (WFT) transform methods [3–6], spatial carrier phase shifting (SCPS) [7–9], regularized phase tracking (RPT) [10–14], continuous wavelet (CWT) [15–20] and S-transform (ST) [21, 22] methods, and Hilbert transform (HT) aided by empirical mode decomposition (EMD) [23–27]. Because of quite extensive literature related to each mentioned technique we have quoted only a few selected papers including original contributions and articles with extensive reference lists. Yet another recent development, influential on the presented paper, is the algorithm of the radial basis functions (thin plate splines in particular) phase interpolation on the fringe pattern skeleton [28]. Two main phase estimation solutions can be distinguished (see, for example, [5, 6, 17, 18]):

Phase estimation: in the WFT, CWT and ST methods the phase is estimated locally from the transform result. This approach provides high estimation accuracy but requires further unwrapping process;
Frequency estimation: the instantaneous frequencies are estimated and subsequently integrated (FT, WFT, CWT, ST). Since the phase distribution obtained in this way is continuous, no further unwrapping process is required.

The phase unwrapping process, especially when processing and analyzing complex fringe patterns with considerable background and modulation (contrast) variations and high noise contamination (e.g., multiplicative speckle noise inherently met in digital speckle pattern interferometry, DSPI, correlation patterns) is not straightforward. The situation becomes even more complex when dealing with closed fringe patterns.

In this paper we present the new algorithm capable of extracting the continuous phase distribution from the noise corrupted interferogram. The method is based on the smoothing spline concept. Splines are piecewise polynomial functions used abundantly for interpolation and curve/surface design [29]. Spline built with polynomials of degree (highest argument power) n is continuous and has first n − 1 continuous derivatives. The smoothing spline [29, 30] gives opportunity to enforce increased smoothness while limiting admissible error, which is of great importance when the noisy data approximation is to be considered. The smoothing spline is not unknown to the optical community, e.g., as a tool for the signal denoising [31, 32]. Overall, it is a well established and very successful concept in signal (and image) processing.

In this work, instead of tackling the task of the noise removal from the intensity fringe image as an image enhancement technique, we address the direct estimation of the phase distribution with the newly introduced algorithm named implicit smoothing spline (ISS). In this way the error in the sought phase distribution can be minimized directly, unlike in the methods aiming at the intensity pattern filtering. The method also offers, under proper conditions (initial phase guess), capability to estimate the continuous phase distribution, without the mod 2π ambiguity or function symmetry ambiguity. Thus, we may avoid potentially difficult phase unwrapping at the cost of providing a reasonable initial approximation to the phase distribution. The common problem of numerous fringe processing algorithms (RPT, CWT, Fourier) in the low fringe frequency regions is not severe for the ISS, since it does not explicitly depend on the concept of the fringe frequency. Finally, it constitutes a new framework for the phase demodulation and with very small number of parameters that need to be specified, the method is not difficult to automatize.

2. About the method

The implicit smoothing spline (ISS) method, reported in this paper, can be regarded as the generalization of the classic smoothing spline algorithm [30]. It allows to define the smoothing spline approximation in the case when the provided noised input data represents some known function of the approximated distribution rather than the approximated distribution itself. In this section we show the formal relation between the smoothing spline and the ISS. We exclusively discuss the cubic spline, as this is by far the most popular variant, well motivated by the cubic spline strain energy minimization property.

2.1. Cubic smoothing spline

The reason behind the significant success of smoothing spline as the data analysis tool is that it offers a simple framework for the balance between smoothness and interpolation error magnitude, regulated with a single parameter, hereafter denoted with p. On the opposite sides of the spectrum of solutions described by the smoothing spline are the direct spline interpolation (no restriction on smoothness, p = 0) and straight line least squares fit (for p → ∞). Formally, the cubic smoothing spline is a spline function of degree n = 3 minimizing the following penalized sum of squares (PSS) functional

P (S) = \sum_{i = 0}^{N} {[Y_{i} - S_{i}]}^{2} + p \int_{Ω} {[S^{″} (x)]}^{2} d x,

where S_i = S(x_i), the spline function S(x) evaluated at the interpolation knots x_i (one knot for each pixel in this case) and Ω represents the interpolation domain, i.e., in 1D it is the interval Ω = [x₀, x_N]. Data value at the i-th interpolation knot is denoted by Y_i. The stunning fact is that there exists an unique solution Ŝ(x), specified with the vector of values Ŝ, minimizing such a PSS functional and that formula for calculating such a spline function can be given explicitly. There is no need for any special numerical procedure aimed at the PSS functional minimization, the solution is found just by the means of the elementary linear algebra. This is possible because the penalizing part of Eq. (1) can be written as a simple bilinear form [30]

\int_{Ω} {[S^{″} (x)]}^{2} d x = S^{T} K S .

K is an N × N matrix of constant coefficients that only depends on the spline degree and the knots distributions (see Appendix A for details in case of the cubic spline). Differentiating the PSS functional we find the system of equations

\frac{\partial P (S)}{\partial S_{i}} = - 2 (Y - S) + 2 p K S .

Equaling derivatives to zero we find the solution Ŝ(x), given by

(I + p K) \hat{S} = Y,

with identity matrix I. The numerical implementations of the smoothing spline typically take advantage of the K matrix sparsity to solve the system of Eqs. (4), which is of particular importance if the large data set is considered.

2.2. Implicit smoothing spline

We aim at exploring the case, in which the data vector Y_i does not correspond directly to the sought distribution S_i = S(x_i), but rather to the known scalar function f of S_i. Therefore we modify the first component (the error) of the PSS functional, leaving the second component (the smoothness) untouched. Following the derivation outlined in the former section, we find

P_{f} (S) = \sum_{i = 0}^{N} {[Y_{i} - f (S_{i})]}^{2} + p \int_{Ω} {[S^{″} (x)]}^{2} d x .

After differentiation and some straightforward algebra we find the following system of algebraic equations

f^{'} (\hat{S}) \cdot [Y - f (\hat{S})] = p K \hat{S},

where we denoted elementwise multiplication by dot and differentiation with the apostrophe. Observe that if f(S) is an identity, i.e., f(S_i) = S_i, Eq. (6) reduces to Eq. (4) and we have a standard smoothing spline problem. We emphasize that this is not the differential equation, since we are looking for the vector Ŝ and the function f (and hence its derivatives) is known explicitly. For normalized fringe pattern, which is of our interest, f(S_i) = cos(S_i) and Eq. (6) becomes

sin (2 \hat{S}) - 2 Y \cdot sin (\hat{S}) = 2 p K \hat{S},

constituting a nonlinear system of equations with the unknown vector Ŝ. This is unlike the case of a regular smoothing spline, where the system to be solved was linear, Eq. (4). However, as long as functions involved in Eq. (7) are continuous and smooth, the gradient-based nonlinear algebraic solvers are expected to work fine. In case of f(S_i) = cos(S_i), the minimizer of the PSS given in Eq. (5) can not be unique even for the simple reason of the cosine function 2π periodicity. This means that some attention has to be given to the choice of the initial condition for the nonlinear solver, which we discuss further.

A simple numerical illustration of the ISS algorithm performance and its comparison to the regular smoothing spline is presented in Fig. 1. The data is in the form of cos(ϕ) + N(x), ϕ = x² function spoiled with the additive white Gaussian noise N(x) of standard deviation σ = 0.5. Figures 1(a) and 1(b) show the phase demodulation approach with the smoothing spline, i.e., denoise the intensity first and then decode phase with direct inversion of the cos(ϕ) function. On the other hand, Figs. 1(c) and 1(d) show the ISS approach where it is phase distribution to be directly approximated by the smoothing spline. Constant function was used as the initial phase guess for the nonlinear solver. While the phase estimation quality obviously depends on the smoothing parameter p, we have repeatedly found the direct phase estimation superior to intensity filtration and subsequent phase decoding.

Fig. 1 Smoothing spline intensity denoising (a) and phase calculation (b) compared with the ISS-estimated phase (d) and the corresponding intensity (c).

Download Full Size | PDF

3. Application for the fringe pattern phase decoding

In order for the ISS algorithm to be applied to the phase decoding of the fringe patterns with f (S_i) = cos(S_i), some details of the method need to be further clarified. In Fig. 2 we show the flowchart of the general algorithm allowing to utilize ISS for the fringe pattern phase decoding. Subsequently, the algorithm is discussed in more details.

Fig. 2 Flowchart of the proposed phase demodulation algorithm based on the ISS.

Download Full Size | PDF

3.1. Nonlinear solver initial condition

Solving the nonlinear system of Eqs. 7 can be problematic without the reasonable initial condition, as distinct parts of the Ŝ(x) may converge to the solution branches mutually shifted by 2kπ, integer k. Initial phase guess must be provided. While the guess does not need to be very accurate and constant value was enough in the example shown in Fig. 1, it should indicate the closest solution branch to the nonlinear solver. In any situation, when the approximate phase distribution is known a priori (e.g., when departure from the ideal object shape is measured), it can be used as the initial guess. In the subsequent demonstrations we address the more complicated case, when nothing is known about the expected phase distribution. Following [28], we use the thin plate spline (TPS) interpolation on the interferogram skeleton set to approximate phase distribution. Such a smooth interpolation built on the ungridded skeleton set is subsequently used as the ISS initial phase guess. This part is indicated with dashed lines on the flowchart.

3.2. Choice of the smoothing parameter

In the ISS algorithm there is a single parameter p, controlling the smoothness of the phase estimate Ŝ(x). Wrong choice of p may result either in preserving the unwanted noise component or oversmoothing the phase estimation. Many methods of selecting the smoothing parameter value were developed for the classic spline interpolation, among which the generalized cross-validation (GCV, see, e.g., [33]) seems to be the most successful. While the basic idea could be directly implemented to ISS, this would result in very much unpractical processing times. One of the reasons of the GCV method success is that effective, non-direct implementation is available. It can only be employed for the classic smoothing spline though, where linear equations are solved. While further investigations can reveal GCV methods (heuristic, most likely) available for the ISS algorithm, we are currently satisfied with the intuitive selection of the smoothness parameter p. Larger values of p provide better estimation quality for the data of lower signal to noise ratio. In all of the numerical tests presented further values of p were kept around p = 0.1.

3.3. Equation or functional?

To find the ISS problem solution we need to solve a nonlinear set of equations (7). On the other hand, we may also attempt to minimize the P_f (S) functional from Eq. (5). Not only these formulations are not strictly equivalent from the mathematical point of view, they also need numerical methods to be employed in a different manners. We investigated both cases, utilizing the powerful nonlinear solver based on the trust-region dogleg method [34, 35] (free Fortran and commercial Matlab implementations available). While for both problems solution should converge to the same vector Ŝ under proper initial conditions, we observed different numerical performance. The method solving Eq. (7) repeatedly found the solution of higher quality, indicating better robustness of this formulation. This observation is consistent with the computational mathematic practice, see chapters IX-X of [36]. Since both approaches yielded similar time performance, we limit further demonstrations to the variant of the ISS implementation in which Eq. (7) is solved.

3.4. Extension to 2D

The examples shown so far were calculated for the single dimensional case. In principle 1D algorithms can be used for the interferometric pattern processing purposes, being applied to each line of pixels separately. Nevertheless, 2D algorithms are preferred for their overall better robustness. Generalizing the presented algorithm to 2D (or, for that matter, arbitrary dimension) is straightforward using the tensor product of the 1D splines in their basis form, as long as the knots are distributed on the regular grid. The details, which are similar for ISS and regular smoothing spline, are thoroughly discussed in the monograph [29], chapter XVII. Nevertheless, we found different approach to be of lower computational complexity and sufficient robustness. We calculate two 1D ISS rasters, for rows and for columns. Then we average so-obtained images and apply standard 2D tensored smoothing spline to the result. In further demonstrations all ISS results refer to such a variant of the algorithm. In the first of the simulations we also show the performance of the method without final 2D smoothing.

3.5. Data normalization

ISS assumes that the normalized interferometric pattern is given as an input. Assuming that the experimentally registered intensity pattern obeys (random noise ignored)

I (x) = b (x) + m (x) cos [ϕ (x)],

with background illumination intensity b(x) and fringe modulation m(x), it is necessary to transform data into form

\hat{I} (x) \approx cos [ϕ (x)] .

Therefore, preprocessing is necessary to remove background illumination influence and the possible contrast variations from the input. However, the noise component does not need to be effectively treated at the preprocessing stage - the filtering properties of the implicit smoothing spline account for it. There are many methods available in the literature to perform the necessary preprocessing when the low quality fringe pattern is analyzed. One example is the fast and data-driven technique based on the fast and adaptive empirical mode decomposition [26, 27], with Matlab implementations available online. If the phase initial guess is based on the fringe skeleton estimation, one can utilize its construction for the pattern normalization, simply transforming the intensity in such a way that its mean value becomes 1 along each top fringe ridge and −1 along each bottom fringe ridge. Note that this very basic procedure also approximately removes the background component.

3.6. Relation to other phase estimation methods

While there exist other algorithms aiming at direct phase estimation, such as the RPT algorithms family, there are fundamental differences between the proposed technique and other methods considered in the literature, namely

phase estimation is typically performed based on the local model, typically of a linear [13] or quadratic [14] polynomial. WFT, CWT and ST can also be regarded as methods employing locally the linear phase model. Spline, on the other hand is a global function represented with a different polynomial between each adjacent knots. Since in our approach every single pixel constitutes a knot, we can represent any function defined on our discrete domain using splines, there is no model limitation whatsoever;
thanks to the smoothing splines algebraic properties, Eq. (2), we globally optimize just the phase values themselves rather than local sets of polynomial coefficients or other model parameters and since we are able to cast the optimization problem in the closed form of the system of algebraic equations depending only on the data values and unknown phase distribution, we may easily apply robust nonlinear algebraic solvers;
global character means stronger dependence on the quality of the initial phase condition than in case of the local approaches. On the other hand this is the reason why the continuous phase distribution is found without the frequency integration;
the smoothness term of the ISS is curvature-based, while methods such as RPT typically enforce smoothness based on the phase variability. This means that if smoothness is pushed to the limit, ISS is likely to find the best-fit constant frequency term while other methods will tend to reduce any phase variability. This implies that the choice of the smoothness parameter p is less critical in the presented approach.

4. Numerical tests

We present processing results for several typical interferometric data inputs. Some of these include regions of very low fringe frequency. All the synthetic data is in 256 × 256 pixels format. The quality of the retrieved phase is evaluated with the root mean square error (RMS), mean absolute value of the error (MAE) and the maximum error value (Peak), all given in radians.

4.1. General comparison between the compared methods

We compare the ISS approach with the wavelet transform method and the skeleton interpolation method of [28]. The former is well-established and successful fringe pattern analysis technique while the latter is somewhat close to the ISS method, which can be regarded as its extension (since in the presented simulations we use the skeleton-interpolated phase as the nonlinear solver initial condition). Below we give short summary of these two methods:

Interpolation on the skeleton is a direct phase estimation based on the set of interpolation knots being the positive and negative fringe ridges, i.e., the skeleton. Thin plate spline (TPS) interpolation is used, since it allows to perform calculations with the ungridded set of knots, results in smooth function and performs well even with low fidelity, disconnected skeleton [28]. The skeleton acts like a set of contour lines of phase distribution, separated by π radians (2π if only positive ridges are taken into account). Thus, the denser are the fringes and higher skeleton evaluation quality, the better method’s performance. Unlike ISS, this approach ignores the points outside of the skeleton. Also, the skeleton is fixed as the set of interpolation knots, while ISS allows, to certain extent, to correct the skeleton calculation errors, as it is only utilized for the nonlinear solver initial guess. Thin plate spline is expensive to evaluate, so this may not be the preferred technique when the skeleton is a very large set, e.g., for the spatial carrier fringe pattern. For some details regarding the skeleton calculation methods, see [2]. In the examples considered further, we assume that fringe ridges are calculated exactly and the adequate integer multiplicity of π is assumed as a ridge value. While this is not the true value at the ridge pixel for the realistic (i.e., discretely sampled) image, it is still a high quality estimation - more than one could expect from the skeleton estimation on the real, noisy data. Note that such an approach renders results independent of the synthetic pattern noise level. We refer to this method further with the TPS acronym. Note that whenever the ISS algorithm utilizes the TPS method output as an initial phase, the final phase error can not become any larger in result of the ISS processing. Hence, to be precise, we rather observe how much more can ISS reduce the phase error than compare ISS to TPS.
Wavelet transform methods constitute both popular and highly acknowledged phase decoding algorithms. In this comparison we utilized 2D complex continuous Morlet wavelet (CWT), that can be regarded as a windowed Fourier transform improved by utilizing additional relation between the frequency and the window width [17]. When it helped the phase calculation quality, we used unusually small value of the envelope wideness parameter (such an effect is not surprising for the low fringe densities). While theoretically this limits the mid-pass filter properties of the CWT, we found overall results to be much better - for typical width of the envelope Morlet CWT could not compete neither with TPS nor ISS in any of the considered examples other than the high frequency carrier image. Big advantage of the CWT is that there is no need for the additional preprocessing such as normalization or skeleton calculation. On the other hand, the resultant phase is mod 2π wrapped. Also, even function ambiguity (cos(ϕ) = cos(−ϕ)) is present, leading to further unwrapping errors. One possible method to deal with this problem is described in [37]. In the error calculations given further we assumed that all wrapped phase ambiguities can be resolved with perfect accuracy. The quality was always worse when we performed the actual phase unwrapping. In the corresponding figures we show the phase decoding results in the wrapped form and error distributions of the perfectly unwrapped data.

4.2. Synthetic interferometric patterns

4.2.1. Test 1: small phase variation

In this example there is no carrier and phase variation is too small to produce fringes, hence we need to simply invert the cosine function in the properly normalized image. The TPS method fails since there are no fringe ridges present. Since no carrier is present, Fourier method can not be applied and the Morlet wavelet application is cumbersome. One very basic approach is to smooth the intensity image and find its arccos function. While there is a plethora of fringe image filtration methods, we only consider the optimal linear filtration, i.e., knowing the noiseless image we choose the Gaussian convolution filter of size that minimizes the error (obviously this is not even doable for the real data, the intention here is to compare ISS performance with some very effective filtration algorithm). The comparison was performed for the phase distribution ϕ₁ shown in Fig. 3(e) for the intensity images with additive white Gaussian noise of standard deviation σ₁ = 0.1 and σ₂ = 0.5. Image of the ϕ₁(x, y) function is confined to the [0, 3] set. Constant phase was utilized as the initial phase guess for the ISS nonlinear solver, since the skeleton interpolation was not feasible. We show phase estimation results for σ₂ in Fig. 3, while phase errors for both noise levels are given in Table 1. While we do not indicate this in the results, it was also possible to perform the phase demodulation with the Morlet 2D CWT, with the strongly narrowed envelope. The error magnitude was similar as in case of the Gaussian filtration, far worse than for the ISS method.

Fig. 3 Synthetic pattern spoiled with noise (a); underlying phase distribution (e); result and residual error of Gaussian filtration and cosine inversion result (b),(f); ISS averaged 1D rasters result and error (c),(g); 2D ISS result and error (d),(h). Errors in radians.

Download Full Size | PDF

Table 1. Errors comparison for the numerical tests 1 and 2

View Table | View all tables in this article

4.2.2. Test 2: continuous phase estimation

Here we consider similar, yet scaled, phase distribution ϕ₂ = 5.25 · ϕ₁. Phase variation is large enough to produce some fringes, see Fig. 4, so the wavelet method is more applicable than in the previous example. TPS skeleton interpolation and ISS are also viable. The Fourier method can not be applied because of the fringe orientation variation. Because of the large fringe curvature in relation to its frequency, envelope width parameter of the CWT was tuned for the higher quality (it is 0.35 of the typically used value). We show results for the noise standard deviation σ₂ = 0.5 in Fig. 4 and give errors for σ₁ = 0.1 and σ₂ = 0.5 in Table 1.

Fig. 4 Synthetic pattern spoiled with noise (a); pattern skeleton used for the interpolation (e); result and error of the TPS method (b),(f); result and residual error of the CWT method (c),(g); ISS result and error (d),(h). Errors in radians.

Download Full Size | PDF

4.2.3. Test 3: fringes with spatial carrier

In this example high frequency vertical carrier fringes distorted by ϕ₃ = 2.625 · ϕ₁ distribution are present. This is the kind of input for which Fourier (classic method described in [3], here referred to as FFT) and CWT methods can be successfully applied. It is rather unusual to address this kind of input with the TPS algorithm, but it can be done. The application of the FFT method was preceded by the linear filtration with Gaussian convolution mask of optimized width, standard deviation equal to 1 pixel and 2 pixels for noise of σ₁ = 0.1 and σ₂ = 0.5, respectively. As previously, figures are shown for the noise of σ₂ = 0.5, see Fig. 5, while the errors for σ₁ = 0.1 and σ₂ = 0.5 are indicated in Table 2.

Fig. 5 Synthetic pattern spoiled with noise (a); ideal noiseless intensity distribution (f); result and error of the FFT method (b),(g); result and residual error of the TPS method (c),(h); result and residual error of the CWT method (d),(i); ISS result and error (e),(j). Errors in radians.

Download Full Size | PDF

Table 2. Errors comparison for the tests 3–5

View Table | View all tables in this article

4.3. DSPI images processing

The multiplicative character of the speckle noise makes DSPI image more challenging for the fringe pattern processing algorithm to decode. Since the mean value of the noise-spoiled intensity is proportional to the sought intensity [38] it is a common practice to analyze speckle pattern with simple low-pass filtration denoising, yet many more sophisticated methods were introduced as well. Here we demonstrate that the ISS algorithm works fine with the speckle fringe images, maintaining all of its benefits. The subtractive speckle fringe pattern may be put in the form

{[Δ I (x)]}^{2} = {[I_{2} (x) - I_{1} (x)]}^{2} \propto {[sin (θ) cos (ϕ / 2)]}^{2} = 0.5 {sin}^{2} (θ) [1 + cos (ϕ)],

where intensities I₂ and I₁ represent speckle patterns registered for two different states of the investigated object, θ is the random speckle phase distribution and ϕ is the sought phase distribution related to the measured physical quantity. We assumed uncorrelated noise values in adjacent pixels (speckles of 1 pixel size), hence the simplified form of Eq. (10). In the ISS framework we may estimate the function of phase f (ϕ) = 1 + cos(ϕ), which, using Eq. (6), leads to the following system of equations

sin (2 \hat{S}) - 2 (Y - 1) \cdot sin (\hat{S}) = 2 p K \hat{S},

Note that if we put Y := Y − 1, we arrive at Eq. (7) again, meaning that the phase estimation in the speckle pattern may simply be performed with the same algorithm as for regular fringe pattern.

4.3.1. Test 4: curved fringes

In this example we consider the same phase distribution distribution as in the test 2, but spoiled with the multiplicative speckle noise, Fig. 6. Corresponding errors for the TPS, CWT and the ISS algorithms are given in Table 2.

Fig. 6 Synthetic DSPI pattern (a); ideal noiseless intensity distribution (e); result and error of the TPS method (b),(f); result and residual error of the CWT method (c),(g); ISS result and error (d),(h). Errors in radians.

Download Full Size | PDF

4.3.2. Test 5: low frequency phase variation

The final synthetic data example involves low frequency speckle pattern distribution, Fig. 7. Errors are given in Table 2.

Fig. 7 Synthetic DSPI pattern (a); ideal noiseless intensity distribution (e); result and error of the TPS method (b),(f); result and residual error of the CWT method (c),(g); ISS result and error (d),(h). Errors in radians.

Download Full Size | PDF

4.4. Processing experimental patterns

Here we consider the experimental DSPI pattern from [39] obtained in nodestructive testing with thermal waves. A temporally modulated IR lamp was used to heat the surface of a flawed plate and the resulting out-of-plane displacements were monitored. There is no background in the data, since it is obtained with the subtraction of two frames. The skeleton is estimated based on the basic binarization and morphological operations. For normalization we only apply the linear transformation based on the estimated skeleton, as indicated in Subsection 3.5. Hence, no nonlinear operation on the image intensities is performed and the ISS algorithm optimizes phase distribution using the image shown in Fig. 8(a). As the phase error distribution is not available for the real data, we decided to present the rewrapped phase distributions, i.e., the cosine function of the calculated phase in order to demonstrate the phase decoding quality. Results are compared in Fig. 8, clearly indicating that the ISS yields smoother phase distribution approximation.

Fig. 8 Experimental DSPI pattern [39] as used by the ISS algorithm (a); enhanced contrast DSPI pattern, shown solely for demonstratory purposes (b); rewrapped CWT-processed pattern (c); rewrapped ISS-processed pattern (d).

Download Full Size | PDF

In Fig. 9 an image obtained in testing a notch including specimen under tensile load is shown. In case of these data the CWT method failed to find a reasonable approximation to the fringe pattern. This is because of the presence of a discontinuity (very much local feature) in the low frequency fringes (large scale changes). There is simply no Morlet wavelet family element that could approximate the intensity in the discontinuity region. Similar problems may appear for any algorithm utilizing locally linear phase model. ISS, on the other hand, performs well in the whole domain. The last experiment, Fig. 10, utilizes the interferometric pattern produced by the spherical mirror tested with Burch’s common-path scatterplate interferometer, [40]. This data is critically difficult for the single-frame algorithm to process since the presence of the doubly unscattered and doubly scattered parasitic terms. We approach it with the basic normalization and background removal mentioned in the Subsection 3.5. We were unable to produce any meaningful result with the CWT method.

Fig. 9 Experimental DSPI pattern as used by the ISS algorithm (a); rewrapped CWT result (b); rewrapped ISS-processed pattern (c).

Download Full Size | PDF

Fig. 10 Experimental interferometic pattern [40], as used by the ISS algorithm (a); rewrapped ISS-processed pattern (b); partially found pattern skeleton used to produce ISS initial phase guess (c).

Download Full Size | PDF

4.5. Remarks on the numerical performance

The ISS algorithm indicates very good performance in all of the considered examples, without doubts yielding quality superior to other discussed algorithms. A single disadvantage of the method is a rather high processing time, in considered cases typically about 3 minutes per image with 3.2GHz CPU personal computer. For the same data CWT or TPS processing typically lasted no more than about 1 minute. While the performance should be significantly improved with more sophisticated implementation, it also depends strongly on the inner parameters of the solver. For instance, changing the maximum number of solver iterations from 100 to 15 resulted in processing time reduced from above 10 minutes to below 3 minutes with accuracy diminished just by about 1 percent. One may attempt to find a reasonable trade-off between the processing speed and quality by adjusting maximum iterations number or other solver parameters. More attention shall be given to this issue in the future. On the other hand one could reduce the number of the spline function knots, which is equal to the number of data points (image pixels) in the discussed algorithm. On the one hand this could dramatically improve the speed of such an implementation. Yet the price would obviously be the loss of the model generality and corresponding resolution and accuracy loss.

5. Conclusions

We have proposed a new data approximation algorithm that can be regarded as the smoothing spline generalization, namely the implicit smoothing spline method. The application for the interferometric fringe pattern analysis was thoroughly discussed and demonstrated. We showed that in the ISS framework, the solution to the phase evaluation problem can be found by solving a set of algebraic equations rather than minimizing the cost functional and found such an approach to be more robust. A variety of numerical examples was given and all of them corroborated the superior quality of the ISS phase decoding results. The largest benefits of the method were expected in case of low density fringes, where most other methods will fail to properly decode small phase variations (see Figs. 7(g)–7(h) and Fig. 9). We found that ISS may outperform other techniques even with the high frequency spatial carrier fringes, Fig. 5(j). Nevertheless, an honest remark should be given, that determining the phase initial guess is expected to be generally more difficult for the higher spatial frequencies of the fringe pattern. Since this paper is the first introduction of the ISS algorithm, we conclude with the brief summary of the further research directions that we see worth following:

smoothing splines were formulated for arbitrary odd polynomial degree [30], in principle one could imagine extension of this algorithm to the quintic (or even higher order) spline interpolation.
weighted errors or spatially varying smoothness parameter (see, e.g., [41]) should work within the ISS framework in the same manner as with the classic smoothing spline. It could be a very useful feature in case when the pattern quality and/or signal to noise ratio vary throughout the domain.
a reasonable method to automatically choose the smoothness parameter p is necessary.
employing (meta)heuristic procedures, such as simulated annealing or genetic algorithms, to solve Eq. (5) would relax the initial phase guess restrictions by enabling the global extremum localization.
more sophisticated ISS implementation will help with the computation time reduction. In particular, since the initial stage of the algorithm works as the set of independent 1D interpolation problems, parallelization should be feasible.

6. Appendix A

Consider cubic spline (n = 3) and the vector of interpolation knots x_i. The differences vector h_i = x_i+1 − x_i. Note that for our needs, when knots correspond simply to subsequent image pixels, h is likely to be constant. Nevertheless, we give a more general equations here. The matrix K, introduced in Eq. (2) and utilized by both smoothing splines and ISS methods, is in the form of

K = Q^{T} G^{- 1} Q

with cubic spline basis Gram matrix of size N − 2 × N − 2

G_{i i} = h_{i} + h_{i + 1}; G_{i - 1, i} = G_{i, i - 1} = \frac{h_{i}}{2}

and matrix Q of size N − 2 × N

Q_{i i} = \frac{1}{h_{i}}; Q_{i, i + 1} = - \frac{h_{i} + h_{i + 1}}{h_{i} h_{i + 1}}; Q_{i, i + 2} = \frac{1}{h_{i + 1}}

While Eq. (4) is correct, from the numerical point of view it is not a preferred way to calculate the smoothing spline. Typically, some equivalent linear equation is solved [29]. For the ISS algorithm implementation, we are currently directly solving numerically Eq. (7).

7. Appendix B

Several definitions of the penalized smoothing spline functional were given in the literature, leading to a confusing definition of the smoothing parameter. All functionals describe the same spectrum of solutions and while finding the relation between different parametrizations is elementary, we believe that readers interested in using smoothing spline or ISS could save some time and avoid confusion using the following table, Table 3, indicating the relation between three popular definitions of the PSS.

{PSS}_{1} (S) = \sum_{i = 1}^{N} {[Y_{i} - S_{i}]}^{2} + p \int_{Ω} {[S^{″} (x)]}^{2} d x

{PSS}_{2} (S) = ρ \sum_{i = 1}^{N} {[Y_{i} - S_{i}]}^{2} + \int_{Ω} {[S^{″} (x)]}^{2} d x

{PSS}_{3} (S) = λ \sum_{i = 1}^{N} {[Y_{i} - S_{i}]}^{2} + (1 - λ) \int_{Ω} {[S^{″} (x)]}^{2} d x

Table 3. Three popular PSS parametrizations

View Table | View all tables in this article

Even more confusion is introduced when matrices given explicitly in the Appendix A differ by a scalar value. Unfortunately, it does happen between different implementations and can not be solved with a simple table.

Acknowledgments

We would like to thank P. Bechler for his comments on the numerical procedures, L. Diet-rich and L. Salbut for access to experimental data and anonymous reviewers for their helpful comments. MW also thanks G. H. Kaufmann and the Institute of Physics in Rosario for their kind hospitality during MW visit to Argentina, where most of the presented research was conducted. This work was supported by the Polish National Science Center grant DEC-2012/07/B/ST7/01392 and the statutory funds of the Faculty of Mechatronics, Warsaw University of Technology. MW also acknowledges the support of the European Union in the framework of European Social Fund through the Warsaw University of Technology Development Programme UE CAS.

References and links

1. D. Robinson and G. Reid, eds., Interferogram Analysis: Digital Fringe Pattern Measurements (Institute of Physics, 1993).

2. D. Malacara, M. Servin, and Z. Malacara, Interferogram Analysis for Optical Testing (CRC, 2005). [CrossRef]

3. M. Takeda, H. Ina, and S. Kobayashi, “Fourier transform methods of fringe-pattern analysis for computer-based topography and interferometry,” J. Opt. Soc. Am. 72(1), 156–160 (1982). [CrossRef]

4. C. Rodier and F. Roddier, “Interferogram analysis using Fourier transform techniques,” Appl. Opt. 26(9), 1668–1673 (1987). [CrossRef]

5. Q. Kemao, “Two-dimensional windowed Fourier transform for fringe pattern analysis: principle, applications and implementations,” Opt. Lasers Eng. 45(2), 304–317 (2007). [CrossRef]

6. S. Fernandez, M. A. Gdeisat, J. Salvi, and D. Burton, “Automatic window size selection in Windowed Fourier Transform for 3D reconstruction using adapted mother wavelets,” Opt. Commun. 284(12), 2797–2807 (2011). [CrossRef]

7. Y. Ichioka and N. Inuiya, “Direct phase detecting system,” Appl. Opt. 11(7), 1507–1514 (1972). [CrossRef] [PubMed]

8. O. Y. Kwon and D. M. Sough, “Multichannel grating phase-shift interferometers,” Proc. SPIE 599, 273–279 (1985). [CrossRef]

9. K. Creath and J. Schmit, “N-point spatial phase-measurement techniques for non-destructive testing,” Opt. Lasers Eng. 24, 365–379 (1996). [CrossRef]

10. M. Servin, J. L. Marroquin, and F. J. Cuevas, “Demodulation of a single interferogram by use of a two-dimensional regularized phase-tracking technique,” Appl. Opt. 36(19), 4540–4548 (1997). [CrossRef] [PubMed]

11. H. Wang and Q. Kemao, “Frequency guided methods for demodulation of a single fringe pattern,” Opt. Express 17(17), 15118–15127 (2009). [CrossRef] [PubMed]

12. C. Tian, Y. Y. Yang, D. Liu, Y. J. Luo, and Y. M. Zhuo, “Demodulation of a single complex fringe interferogram with a path-independent regularized phase-tracking technique,” Appl. Opt. 49(2), 170–179 (2010). [CrossRef] [PubMed]

13. L. Kai and Q. Kemao, “A generalized regularized phase tracker for demodulation of a single fringe pattern,” Opt. Express 20(11), 12579–12592 (2012). [CrossRef] [PubMed]

14. L. Kai and Q. Kemao, “Improved generalized regularized phase tracker for demodulation of a single fringe pattern,” Opt. Express 21(20), 24385–24397 (2013). [CrossRef] [PubMed]

15. L. Watkins, S. Tan, and T. Barnes, “Determination of interferometer phase distributions by use of wavelets,” Opt. Lett. 24(13), 905–907 (1999). [CrossRef]

16. Z. Wang and H. Ma, “Advanced continuous wavelet transform algorithm for digital interferogram analysis and processing,” Opt. Eng. 45(4), 045601 (2006). [CrossRef]

17. M. A. Gdeisat, D. R. Burton, and D. R. Lalor, “Spatial carrier fringe pattern demodulation by use of a two-dimensional continuous wavelet transform,” Appl. Opt. 45(34), 8722–8732 (2006). [CrossRef] [PubMed]

18. L. R. Watkins, “Review of fringe pattern phase recovery using 1-D and 2-D continuous wavelet transforms,” Opt. Lasers Eng. 50(8), 1015–1022 (2012). [CrossRef]

19. K. Pokorski and K. Patorski, “Processing and phase analysis of fringe patterns with contrast reversals,” Opt. Express 21(19), 22596–22614 (2013). [CrossRef] [PubMed]

20. P. Etchepareborda, A. L. Vadnjal, A. Federico, and G. H. Kaufmann, “Phase-recovery improvement using analytic wavelet transform analysis of a noisy interferogram cepstrum,” Opt. Lett. 37(18), 3843–3845 (2012). [CrossRef] [PubMed]

21. A. Dursun, Z. Sarac, H. S. Topkara, S. Ozder, and F. N. Ecevit, “Phase recovery from interference fringes by using S-transform,” Measurement 41(4), 403–411 (2008). [CrossRef]

22. M. Zhong, W. Chen, T. Wang, and X. Su, “Application of two-dimensional S-Transform in fringe pattern analysis,” Opt. Lasers Eng. 51(10), 1138–1142 (2013). [CrossRef]

23. N. E. Huang, Z. Shen, S. R. Long, M. C. Wu, H. H. Shih, Q. Zheng, N.-C. Yen, C. C. Tung, and H. H. Liu, “The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis,” Proc. R. Soc. London A 454, 903–995 (1998). [CrossRef]

24. N. E. Huang and Z. H. Wu, “A review on Hilbert-Huang transform method and its applications to geophysical studies,” Rev. Geophys. 46(2), 1–23 (2008). [CrossRef]

25. Y. Lei, J. Lin, Z. He, and M. Zuo, “A review on empirical mode decomposition in fault diagnosis of rotating machinery,” Mech. Syst. Signal Process. 35(1–2), 108–126 (2013). [CrossRef]

26. M. Trusiak, K. Patorski, and M. Wielgus, “Adaptive enhancement of optical fringe patterns by selective reconstruction using FABEMD algorithm and Hilbert spiral transform,” Opt. Express 20(21), 23463–23479 (2012). [CrossRef] [PubMed]

27. M. Trusiak, M. Wielgus, and K. Patorski, “Advanced processing of optical fringe patterns by automated selective reconstruction and enhanced fast empirical mode decomposition,” Opt. Lasers Eng. 52(1), 230–240 (2014). [CrossRef]

28. G. Wang, Y. J. Li, and H. Ch. 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]

29. C. de Boor, A Practical Guide to Splines (Springer, 1994).

30. C. Reinsch, “Smoothing by spline functions,” Numer. Math. 10, 177–183 (1967). [CrossRef]

31. A. Federico and G. H. Kaufmann, “Local denoising of digital speckle pattern interferometry fringes by multiplicative correlation and weighted smoothing splines,” Appl. Opt. 44(14), 2728–2735 (2005). [CrossRef] [PubMed]

32. M. J. Peyrovian and A. A. Sawchuk, “Restoration of noisy blurred images by a smoothing spline filter,” Appl. Opt. 16(12), 3147–3153 (1977). [CrossRef] [PubMed]

33. P. Craven and G. Wahba, “Smoothing noisy data with spline functions estimating the correct degree of smoothing by the method of generalized cross validation,” Numer. Math. 31, 377–403 (1979). [CrossRef]

34. M. J. D. Powell, “A hybrid method for nonlinear equations,” in Numerical Methods for Nonlinear Algebraic Equations, P. Rabinowitz, ed. (Gordon and Breach, 1970), pp. 87–114.

35. M. J. D. Powell, “A Fortran subroutine for solving systems of nonlinear algebraic equations,” in Numerical Methods for Nonlinear Algebraic Equations, P. Rabinowitz, ed. (Gordon and Breach, 1970), pp. 115–161.

36. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes: The Art of Scientific Computing, 2 (Cambridge University, 1992).

37. J. Ma, Z. Wang, B. Pan, T. Hoang, M. Vo, and L. Luu, “Two-dimensional continuous wavelet transform for phase determination of complex interferograms,” Appl. Opt. 50(16), 2425–2430 (2011). [CrossRef] [PubMed]

38. J. W. Goodman, Speckle Phenomena in Optics: Theory and Applications (Roberts and Company, 2007).

39. G. H. Kaufmann, “Nondestructive testing with thermal waves using phase shifted temporal speckle pattern interferometry,” Opt. Eng. 42(7), 2010–2014 (2003). [CrossRef]

40. K. Patorski and L. Salbut, “Simple polarization phase-stepping scatterplate interferometry,” Opt. Eng. 43(2), 393–397 (2004). [CrossRef]

41. C. de Boor, “Calculation of the smoothing spline with weighted roughness measure,” this paper can be downloaded at http://www.cs.wisc.edu.

	Test 1						Test 2
σ	0.1			0.5			-	0.1		0.5
	Simple	ISS1D	ISS2D	Simple	ISS1D	ISS2D	TPS	CWT	ISS	CWT	ISS
RMS	0.178	0.076	0.056	0.245	0.182	0.136	0.327	0.222	0.124	0.224	0.211
MAE	0.087	0.054	0.043	0.166	0.139	0.115	0.196	0.151	0.083	0.161	0.131
Peak	0.625	0.396	0.198	1.273	0.853	0.436	1.167	1.778	0.214	1.849	1.299

	Test 3							Test 4			Test 5
σ	-	0.1			0.5			-			-
	TPS	FFT	CWT	ISS	FFT	CWT	ISS	TPS	CWT	ISS	TPS	CWT	ISS
RMS	0.282	0.172	0.092	0.057	0.199	0.104	0.075	0.327	0.287	0.253	0.140	0.202	0.083
MAE	0.185	0.162	0.061	0.035	0.173	0.067	0.053	0.196	0.212	0.173	0.089	0.134	0.050
Peak	1.614	0.462	0.496	0.389	0.661	0.563	0.543	1.167	1.570	1.067	0.625	0.974	0.543

-	PSS₁(p) = ...	PSS₂(ρ) = ...	PSS₃(λ) = ...

... = PSS₁(?)	p	$\frac{1}{ρ}$	$\frac{1 - λ}{λ}$
... = PSS₂(?)	$\frac{1}{p}$	ρ	$\frac{λ}{1 - λ}$
... = PSS₃(?)	$\frac{1}{1 + p}$	$\frac{ρ}{1 + ρ}$	λ

	Test 1						Test 2
σ	0.1			0.5			-	0.1		0.5
	Simple	ISS1D	ISS2D	Simple	ISS1D	ISS2D	TPS	CWT	ISS	CWT	ISS
RMS	0.178	0.076	0.056	0.245	0.182	0.136	0.327	0.222	0.124	0.224	0.211
MAE	0.087	0.054	0.043	0.166	0.139	0.115	0.196	0.151	0.083	0.161	0.131
Peak	0.625	0.396	0.198	1.273	0.853	0.436	1.167	1.778	0.214	1.849	1.299

	Test 3							Test 4			Test 5
σ	-	0.1			0.5			-			-
	TPS	FFT	CWT	ISS	FFT	CWT	ISS	TPS	CWT	ISS	TPS	CWT	ISS
RMS	0.282	0.172	0.092	0.057	0.199	0.104	0.075	0.327	0.287	0.253	0.140	0.202	0.083
MAE	0.185	0.162	0.061	0.035	0.173	0.067	0.053	0.196	0.212	0.173	0.089	0.134	0.050
Peak	1.614	0.462	0.496	0.389	0.661	0.563	0.543	1.167	1.570	1.067	0.625	0.974	0.543

Continuous phase estimation from noisy fringe patterns based on the implicit smoothing splines

Abstract

1. Introduction

2. About the method

2.1. Cubic smoothing spline

2.2. Implicit smoothing spline

3. Application for the fringe pattern phase decoding

3.1. Nonlinear solver initial condition

3.2. Choice of the smoothing parameter

3.3. Equation or functional?

3.4. Extension to 2D

3.5. Data normalization

3.6. Relation to other phase estimation methods

4. Numerical tests

4.1. General comparison between the compared methods

4.2. Synthetic interferometric patterns

4.2.1. Test 1: small phase variation

4.2.2. Test 2: continuous phase estimation

4.2.3. Test 3: fringes with spatial carrier

4.3. DSPI images processing

4.3.1. Test 4: curved fringes

4.3.2. Test 5: low frequency phase variation

4.4. Processing experimental patterns

4.5. Remarks on the numerical performance

5. Conclusions

6. Appendix A

7. Appendix B

Acknowledgments

References and links

Cited By

Figures (10)

Tables (3)

Equations (17)

Optics Express