OSA's Digital Library

Optics Express

Optics Express

  • Editor: Andrew M. Weiner
  • Vol. 22, Iss. 9 — May. 5, 2014
  • pp: 10775–10791
« Show journal navigation

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

Maciek Wielgus, Krzysztof Patorski, Pablo Etchepareborda, and Alejandro Federico  »View Author Affiliations


Optics Express, Vol. 22, Issue 9, pp. 10775-10791 (2014)
http://dx.doi.org/10.1364/OE.22.010775


View Full Text Article

Acrobat PDF (2690 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

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

1. Introduction

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

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

]. Spline built with polynomials of degree (highest argument power) n is continuous and has first n − 1 continuous derivatives. The smoothing spline [29

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

, 30

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

] 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

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

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]

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

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)=i=0N[YiSi]2+pΩ[S(x)]2dx,
(1)
where Si = S(xi), the spline function S(x) evaluated at the interpolation knots xi (one knot for each pixel in this case) and Ω represents the interpolation domain, i.e., in 1D it is the interval Ω = [x0, xN]. Data value at the i-th interpolation knot is denoted by Yi. 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

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

]
Ω[S(x)]2dx=STKS.
(2)
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)Si=2(YS)+2pKS.
(3)
Equaling derivatives to zero we find the solution Ŝ(x), given by
(I+pK)S^=Y,
(4)
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 Yi does not correspond directly to the sought distribution Si = S(xi), but rather to the known scalar function f of Si. 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
Pf(S)=i=0N[Yif(Si)]2+pΩ[S(x)]2dx.
(5)
After differentiation and some straightforward algebra we find the following system of algebraic equations
f(S^)[Yf(S^)]=pKS^,
(6)
where we denoted elementwise multiplication by dot and differentiation with the apostrophe. Observe that if f(S) is an identity, i.e., f(Si) = Si, 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(Si) = cos(Si) and Eq. (6) becomes
sin(2S^)2Ysin(S^)=2pKS^,
(7)
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(Si) = cos(Si), 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), ϕ = x2 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).

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 (Si) = cos(Si), 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.

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 2, 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

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]

], 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

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]

]) 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 Pf (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

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

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.

] (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

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

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

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

], 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)],
(8)
with background illumination intensity b(x) and fringe modulation m(x), it is necessary to transform data into form
I^(x)cos[ϕ(x)].
(9)
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

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

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]

], 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

    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]

    ] or quadratic [14

    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]

    ] 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

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

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.

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

table-icon
View This Table
| View All Tables

4.2.2. Test 2: continuous phase estimation

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.

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

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]

], 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, respectively. As previously, figures are shown for the noise of σ2 = 0.5, see Fig. 5, while the errors for σ1 = 0.1 and σ2 = 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.

Table 2. Errors comparison for the tests 3–5

table-icon
View This Table
| View All Tables

4.3. DSPI images processing

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.

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.

4.4. Processing experimental patterns

Here we consider the experimental DSPI pattern from [39

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

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

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

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

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

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

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

    ], 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

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

    ]) 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 xi. The differences vector hi = xi+1xi. 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=QTG1Q
(12)
with cubic spline basis Gram matrix of size N − 2 × N − 2
Gii=hi+hi+1;Gi1,i=Gi,i1=hi2
(13)
and matrix Q of size N − 2 × N
Qii=1hi;Qi,i+1=hi+hi+1hihi+1;Qi,i+2=1hi+1
(14)
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

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

]. For the ISS algorithm implementation, we are currently directly solving numerically Eq. (7).

7. Appendix B

Table 3. Three popular PSS parametrizations

table-icon
View This Table
| View All Tables

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.

OCIS Codes
(100.5070) Image processing : Phase retrieval
(120.2650) Instrumentation, measurement, and metrology : Fringe analysis
(120.3940) Instrumentation, measurement, and metrology : Metrology
(120.5050) Instrumentation, measurement, and metrology : Phase measurement
(120.6160) Instrumentation, measurement, and metrology : Speckle interferometry

ToC Category:
Instrumentation, Measurement, and Metrology

History
Original Manuscript: January 23, 2014
Revised Manuscript: April 8, 2014
Manuscript Accepted: April 10, 2014
Published: April 28, 2014

Citation
Maciek Wielgus, Krzysztof Patorski, Pablo Etchepareborda, and Alejandro Federico, "Continuous phase estimation from noisy fringe patterns based on the implicit smoothing splines," Opt. Express 22, 10775-10791 (2014)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-22-9-10775


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. D. Robinson, G. Reid, eds., Interferogram Analysis: Digital Fringe Pattern Measurements (Institute of Physics, 1993).
  2. D. Malacara, M. Servin, Z. Malacara, Interferogram Analysis for Optical Testing (CRC, 2005). [CrossRef]
  3. M. Takeda, H. Ina, 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, 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, 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, N. Inuiya, “Direct phase detecting system,” Appl. Opt. 11(7), 1507–1514 (1972). [CrossRef] [PubMed]
  8. O. Y. Kwon, D. M. Sough, “Multichannel grating phase-shift interferometers,” Proc. SPIE 599, 273–279 (1985). [CrossRef]
  9. K. Creath, 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, 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, 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, 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, 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, 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, T. Barnes, “Determination of interferometer phase distributions by use of wavelets,” Opt. Lett. 24(13), 905–907 (1999). [CrossRef]
  16. Z. Wang, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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 .

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited