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

We introduce the algorithm for the direct phase estimation from the single noisy interferometric pattern. The method, named implicit smoothing spline (ISS), can be regarded as a formal generalization of the smoothing spline interpolation for the case when the interpolated data is given implicitly. We derive the necessary equations, discuss the properties of the method and address its application for the direct estimation of the continuous phase in both classical interferometry and digital speckle pattern interferometry (DSPI). The numerical illustrations of the algorithm performance are provided to corroborate the high quality of the results. © 2014 Optical Society of America OCIS codes: (120.2650) Fringe analysis; (100.5070) Phase retrieval; (120.5050) Phase measurement; (120.6160) Speckle interferometry; (120.3940) Metrology. 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). 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). 4. C. Rodier and F. Roddier, “Interferogram analysis using Fourier transform techniques,” Appl. Opt. 26(9), 1668– 1673 (1987). 5. Q. Kemao, “Two-dimensional windowed Fourier transform for fringe pattern analysis: principle, applications and implementations,” Opt. Lasers Eng. 45(2), 304–317 (2007). 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). 7. Y. Ichioka and N. Inuiya, “Direct phase detecting system,” Appl. Opt. 11(7), 1507–1514 (1972). 8. O. Y. Kwon and D. M. Sough, “Multichannel grating phase-shift interferometers,” Proc. SPIE 599, 273–279 (1985). 9. K. Creath and J. Schmit, “N-point spatial phase-measurement techniques for non-destructive testing,” Opt. Lasers Eng. 24(5–6), 365–379 (1996). 10. M. Servin, J. L. Marroquin, and F. J. Cuevas, “Demodulation of a single interferogram by use of a twodimensional regularized phase-tracking technique,” Appl. Opt. 36(19), 4540–4548 (1997). 11. H. Wang and Q. Kemao, “Frequency guided methods for demodulation of a single fringe pattern,” Opt. Express 17(17), 15118–15127 (2009). 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). 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). #205384 $15.00 USD Received 23 Jan 2014; revised 8 Apr 2014; accepted 10 Apr 2014; published 28 Apr 2014 (C) 2014 OSA 5 May 2014 | Vol. 22, No. 9 | DOI:10.1364/OE.22.010775 | OPTICS EXPRESS 10775 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). 15. L. Watkins, S. Tan, and T. Barnes, “Determination of interferometer phase distributions by use of wavelets,” Opt. Lett. 24(13), 905–907 (1999). 16. Z. Wang and H. Ma, “Advanced continuous wavelet transform algorithm for digital interferogram analysis and processing,” Opt. Eng. 45(4), 045601 (2006). 17. M. A. Gdeisat, D. R. Burton, and D. R. Lalor, “Spatial carrier fringe pattern demodulation by use of a twodimensional continuous wavelet transform,” Appl. Opt. 45(34), 8722–8732 (2006). 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). 19. K. Pokorski and K. Patorski, “Processing and phase analysis of fringe patterns with contrast reversals,” Opt. Express 21(19), 22596–22614 (2013). 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). 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). 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). 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). 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). 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). 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). 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). 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). 29. C. de Boor, A Practical Guide to Splines (Springer, 1994). 30. C. Reinsch, “Smoothing by spline functions,” Numer. Math. 10, 177–183 (1967). 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). 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). 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). 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, 2nd ed. (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). 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). 40. K. Patorski and L. Salbut, “Simple polarization phase-stepping scatterplate interferometry,” Opt. Eng. 43(2), 393–397 (2004). 41. C. de Boor, “Calculation of the smoothing spline with weighted roughness measure,” this paper can be downloaded at http://www.cs.wisc.edu.


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 pow-erful 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][4][5][6], spatial carrier phase shifting (SCPS) [7][8][9], regularized phase tracking (RPT) [10][11][12][13][14], continuous wavelet (CWT) [15-20] and S-transform (ST) [21,22] methods, and Hilbert transform (HT) aided by empirical mode decomposition (EMD) [23][24][25][26][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.

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.

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 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 0 , 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] 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 ∂ P(S) Equaling derivatives to zero we find the solutionŜ(x), given by 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.

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 After differentiation and some straightforward algebra we find the following system of algebraic equations 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 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 2 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.

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.

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.

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 crossvalidation (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.

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.

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.

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) with background illumination intensity b(x) and fringe modulation m(x), it is necessary to transform data into formÎ 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.

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.

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.

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: 1. 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.
2. 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.

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 φ 1 shown in Fig. 3(e) for the intensity images with additive white Gaussian noise of standard deviation σ 1 = 0.1 and σ 2 = 0.5. Image of the φ 1 (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 σ 2 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.

Test 2: continuous phase estimation
Here we consider similar, yet scaled, phase distribution φ 2 = 5.25 · φ 1 . 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 σ 2 = 0.5 in Fig. 4 and give errors for σ 1 = 0.1 and σ 2 = 0.5 in Table 1.

Test 3: fringes with spatial carrier
In this example high frequency vertical carrier fringes distorted by φ 3 = 2.625 · φ 1 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 σ 1 = 0.1 and σ 2 = 0.5,    Table 2.

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 where intensities I 2 and I 1 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 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.

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.

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.

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. 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.

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.

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). Nev-ertheless, 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.

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 and matrix Q of size N − 2 × N 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).

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.  λ 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.