Hybrid single shot algorithm for accurate phase demodulation of complex fringe patterns is proposed. It employs empirical mode decomposition based adaptive fringe pattern enhancement (i.e., denoising, background removal and amplitude normalization) and subsequent boosted phase demodulation using 2D Hilbert spiral transform aided by the Principal Component Analysis method for novel, correct and accurate local fringe direction map calculation. Robustness to fringe pattern significant noise, uneven background and amplitude modulation as well as local fringe period and shape variations is corroborated by numerical simulations and experiments. Proposed automatic, adaptive, fast and comprehensive fringe analysis solution compares favorably with other previously reported techniques.
© 2016 Optical Society of America
Optical metrology provides noninvasive, full field, easily automatable and accurate tools for engineering/biological object characterization and physical phenomena examination. Applying interferometry (two-beam, speckle, grating, time-averaged etc.), digital holography, optical tomography, structured illumination or moiré techniques one encounters a measurand encoded in phase distribution of a fringe pattern, i.e., fringe local shape, period and position. Measurand information is obtained performing fringe pattern phase demodulation. Multiple-frame methods are recognized as the most accurate ones [1–4 ]. Several fringe patterns describing the same state of measured object have to be acquired under stable conditions in precise manner. This limits the use of multiple-frame techniques to rather slowly varying objects under controlled conditions. Single-frame techniques are very attractive enabling efficient out-of-laboratory tests and accurate studies of transient events in hostile environment. Among dynamically developed single frame automatic fringe pattern analysis (AFPA) techniques we have: Fourier transform [5–7 ], windowed Fourier transform [8–10 ], continuous wavelet transform [11–13 ], regularized phase tracking [14–16 ], spatial carrier phase shifting [17–22 ], statistical image decomposition [23,24 ], empirical mode decomposition based methods [25–31 ] aided by Hilbert spiral transform [32,33 ] and monogenic (analytic/quadrature) signal based approach . The most versatile and efficient solutions have very high computational load and long processing time [13,16,24 ]. Practical challenges include appropriate carrier fringe applications [5–13,17–22 ], error propagation [5–7,14–22 ], susceptibility to noise, fringe modulation/background/period local significant variations [5–24 ], and closed fringes analysis [5–24 ]. Particular case of closed fringes is encountered very often where there is no setup solution for carrier fringe implementation or the phenomena under study are too fast; or a carrier is not introduced for the sake of accuracy. Additionally, in general, fringes with significant orientation and period local variations (considerable horizontal or vertical phase gradient), are especially interesting - classical approaches with the Fourier/Hilbert/wavelet transforms (in general quadrature transforms or analytic signal constructing methods) fail to demodulate phase due to very complicated, overlapping fringe pattern spectrum and sign ambiguity. To enable accurate phase distribution retrieval from a single fringe pattern the information about local fringe direction map is needed. Fringe direction map guides the phase demodulation process.
Fringe orientation is one of the basic characteristics describing a fringe pattern. Following already well-established approach  we define a fringe orientation map as a modulo π indicator. It stores the information about the azimuth of vector locally normal to fringes (e.g., vertical). The fringe direction map is a modulo 2π indicator. It stores the information about an azimuth of vector locally normal to fringes as well as its direction (e.g., up or down). In other words fringe orientation would have the same value at two pixels situated symmetrically to the center of a closed fringe pattern; fringe direction would have two different values (+/− pi). The fringe orientation estimation is cumbersome because we have intensity pattern with phase in the argument of cosine function (and orientation by definition is simply arctangent of the orthogonal spatial derivatives of phase function). The fringe direction map needs even more attention, i.e., advanced unwrapping of the orientation map (from not so informative modulo π to very informative modulo 2π indicator) for spurious angle jumps neglection.
The fringe orientation/direction map is essential in various fringe processing and analysis tasks, i.e., fringe filtering (denoising) [36–55 ], boundary padding [52,55 ], skeletoning (contouring/following/tracking) [38,40,43–45,48,49,51,52,56–58 ], local spatial frequency (fringe period) estimation [41,46,59,60 ] and phase demodulation [39,41,43,44,48,50,58–76 ]. Each of the mentioned operations is guided by a local fringe orientation/direction to enable or enhance the process. In general one can identify several groups of techniques for the fringe orientation estimation (using single fringe pattern): gradient methods [51,57,74,77 ], plane fit methods [42,43 ], combined gradient and plane-fit , spin filters [36,37,39,40,44,45 ] and binary sign-maps [38,40 ], regularized methods [41,53,61–68 ], 2D energy operators , accumulate differences , Fourier transform , Windowed Fourier Transform , Principal Component Analysis based methods [58,65,73 ] and two frame methods, e.g., optical flow .
Fringe orientation map (modulo π) needs to undergo an unwrapping operation to estimate the fringe direction map (modulo 2π) of interest. This operation complicates the direction map determination and lengthens the processing time. In principle it is similar to phase unwrapping issue . Several techniques were employed for this purpose. Quadratic energy function minimization  and regularized phase tracking  are computationally costly due to optimization process. Fringe orientation unwrapping by classifying the gradient vector trends  needs denoising to obtain reliable gradient vectors. Another scheme utilizes quality-guided orientation unwrapping  similar to quality guided phase unwrapping. Advanced and time-consuming Windowed Fourier Transform based processing is required, however.
We propose a novel, comprehensive, all-in-one algorithm for efficient fringe pattern analysis of arbitrary shape and quality. It combines a novel PCA based robust direction estimator with the Hilbert-Huang Transform based fringe analysis technique. Remarkable accuracy of phase demodulation and robustness to noise, background and amplitude variations as well as fringe period and shape changes is obtained taking advantage of local, windowed and global character of the proposed novel hybrid HHT-PCA technique. It is a local scheme due to the pixel-by-pixel automatic selective reconstruction method  employed for very efficient fringe pre-filtering (first step of our approach). It is a windowed scheme due to sliding-window based PCA orientation estimator with the novel segmented-region-based method for robust direction map calculation by area matching (second step of our approach). It is a global scheme due to Larkin's 2D spiral Hilbert transform based direction-aided accurate phase demodulation method (final, third step of our approach). Proposed HHT-PCA scheme is able to accurately perform every fringe processing (denoising, detrending) and analysis task (phase and amplitude demodulation) with robustness to fringe local shape (orientation, spacing) and quality deteriorations (uneven background/modulation and noise influence).
The paper is constructed as follows. Section 2 describes the proposed hybrid HHT-PCA methodology with the emphasis put on the novel robust PCA direction estimator and modified adaptive HHT filtering/demodulation. Such a profitable and original combination of those two novel image processing techniques is reported for the first time. Section 3 contains results of synthetic and experimental evaluation of the proposed scheme employed for the single-shot fringe pattern phase demodulation. The proposed HHT-PCA scheme is compared with other recently reported all-in-one fringe analysis techniques. Comprehensive error analysis is provided. Both quantitative and qualitative evaluation results are discussed in Section 5. Section 6 concludes the paper.
2. Algorithm description
2.1 Principal component analysis based conditioned region matching algorithm for fringe direction map determination
Proposed fringe direction estimator is based on the Principal Component Analysis method and conditioned region matching algorithm. Its working principle will be presented utilizing a simulated fringe pattern, Fig. 1(a) . It has linear background (varying from 0.5 to 2), Gaussian amplitude modulation (varying radially from 1 to 0.2) and multiplicative noise (0.1 variance). First step of the proposed algorithm is the fringe orientation map determination using the results of PCA employed for fringe pattern. We calculate the image gradients in x and y directions using convolution with partial derivative Gaussian kernels . This approach for gradient calculation is robust to fringe pattern noise. Next, using the gradient information, for each pixel of fringe pattern the covariance matrix is calculated
Subsequently employing arc tangent function for orthogonal (x-axis and y-axis) eigenvector components a map of orientation angles (with respect to x-axis) is calculated. To ensure correct processing we assume that eigenvector x-axis component is negative for vertical fringes and y-axis component is positive for horizontal fringes (it ensures efficient orientation vortex detection on later stage of the algorithm). Resulting fringe angle map, presented in Fig. 1(b), contains two types of continuous regions - the first one stores angles in the [1/4π, π] range while the second one in the [-π, −3/4π] range - covering in total a range of π radians. Therefore we have the fringe orientation map which needs to be transformed into a fringe direction map ranging from 0 to 2π radians. A region matching algorithm will be employed for the direction map generation.
Firstly continuous angle regions should be separated, and subsequently depending on the angle values at the region borders the orientation value is corrected. Preliminary continuous region separation is performed in the segmentation process with continuity criterion in terms of constant sign of orientation angle, Fig. 1(c). Obtained regions are not well separated, however. Links are created in locations of the zero phase gradient values.
Next stage of the algorithm consists in efficient zero-phase-gradient mask designing for complete region separation in the orientation angle map. Using preliminary orientation map, Fig. 1(b), we are searching for regions with absolute value of orientation angle larger than 7/8π for horizontal direction (vertical fringes, Fig. 2(a) ) and regions with orientation angle in (3/8π, 5/8π) range for vertical direction (horizontal fringes, Fig. 2(b)). Next two binary images (representing two orthogonal directions) undergo morphological dilation (two times each) and their logical product is computed. Resulting binary image is processed again with dilation (two times) and presented in Fig. 2(c). White regions indicate zero-phase-gradient suspects.
In case of severely defected fringe patterns (including noise, background and contrast variations) in the neighborhood of the zero-phase-gradient spot two locations of undetermined angle (we will call them orientation vortices) are detected, instead of one (specifically located in a zero-phase-gradient spot). Distance between those pixels defines the appropriate mask size for this spot. Orientation vortices are located inside the preliminary masks designed in previous step, Fig. 2(c). Their accurate coordinates can be derived according to the distribution of the relative difference of two eigenvalues of covariance matrix, Fig. 1(d). For each region in Fig. 2(c) its minimum bounding box is calculated. Next coordinates of five starting pixels are computed - 4 pixels at the rectangle corners and fifth one at the diagonals crossing point, Fig. 3(a) . In a 8-neighborhood of each starting pixel we search for minimum value of relative eigenvalue difference, Fig. 3(b), taking into account that this minimum should be located within the zero-phase-gradient suspect region (white pixels in Fig. 3(a)). Starting from each of five pixels we iteratively detect one or more local minima of relative eigenvalue difference. Each minimum is located in a neighborhood of the undetermined angle spot - orientation vortex, Fig. 3(c). The flow-chart of minima detection is depicted in Fig. 4 .
Not every detected region in Fig. 2(c) encompasses a zero-phase-gradient spot. Encountered significant variation of orientation values could be seen as well as the result of edge effects or noise (which is treated as orientation variations; noise “simulates” dense fringes in areas of very low spatial frequency). Characteristic feature of zero-phase-gradient spot in the preliminary orientation map (Fig. 1(b)) is full variability of orientation angles in [0, π) range. In case of minimum (one or more minima) detection in regions determined in Fig. 2(c) we need to design a rectangle (like in Fig. 5 ) to encapsulate it (them) and check how orientation angle values change along its sides. If all orientation values (within assumed subregions of [0, π)) can be encountered then it is an indicator of the zero-phase-gradient spot. Otherwise selected region (i.e., four areas in image corners, Fig. 2(c)) should be ignored (caused by edge effects or noise). Please note that decrease in the signal-to-noise ratio combined with the fringe spatial frequency decrease results in augmented number of regions falsely suspected for encapsulating zero-phase-gradient spots.
Zero-phase-gradient spot detection imposes a necessity to create a mask upon its location. Two cases are to be considered. In the first one rectangle encompasses two regions with negative orientation angle values. They are situated symmetrically with respect to the theoretical undetermined angle value - orientation vortex (Fig. 5(b) and 5(c)). Generated mask should cover all detected minima with mask size minimized. Typically, however, two regions of negative orientation angles are significantly separated hence two locations (instead of one) in a neighborhood of the zero-phase-gradient spot are detected (after operations resulting in Fig. 2(c)). In this case two rectangles will be spanned to encapsulate both negative angle regions (one such a rectangle is shown in Fig. 5(a)). Having all of the spots detected we need to connect the ones located closest to each other (with full variability of orientation angle values and exactly one region of negative orientation angle values). Beside the proximity criterion we should take into account the direction of the orientation angle variation in both locations (clockwise or counter-clockwise) - it should be the same in both detected locations in a neighborhood of the zero-phase-gradient spot. Having all necessary locations connected the masks can be designed (with size minimization criterion in mind).
Fringe orientation map, Fig. 1(b), after determining continuous angle regions with masks is presented in Fig. 6(a) . Next step consists in matching the regions in the orientation map. Values in the largest region (largest number of pixels) are considered as correct. Starting with this region other regions are enclosed starting with the biggest one left (preventing error propagation - the smaller the region the riskier its enclosing is). Orientation angle values in the region to be enclosed are determined by comparing values on the border between two regions under consideration (one with fixed angle values and one to be connected). For every couple of neighboring pixels (4-neighbours) absolute value of the orientation angle difference (θdiff) is calculated as
The result within the (1/2π, 3/2π) range corresponds to negative scalar product of eigenvectors (more computational load) whereas the positive one corresponds to the result in the range (0, 1/2π] [3/2π, 2π). If the second case is more numerous (more positive than negative scalar products) then orientation angles in a region to be enclosed are not to be changed. Otherwise orientation angles should be changed prior to region connection. For negative angle values we add π and from the positive ones π should be subtracted - the conditioned region matching operation transfers fringe orientation to fringe direction.
The last part of the region matching algorithm consists in the direction value determination in pixels previously masked out. For each mask we start by selecting a pixel with the largest relative eigenvalue difference (Fig. 1(d)) located in a 8-neighborhood of a pixel with determined direction angle value. One neighbor (from 8) with fixed direction angle and the largest relative eigenvalues difference is selected. Then for two selected pixels (analyzed one - with undetermined direction angle value; reference one - with determined direction angle value) the absolute value of orientation angle difference is computed, Eq. (3). Value in the analyzed pixel is corrected with the same approach as at the region connection stage. In this way pixels with maximum risk (zero-phase-gradients) are processed last (algorithm proceeds around and towards the orientation vortices). Repeating this scheme for every pixel in every mask we end up with continuous fringe direction map - presented in Fig. 6(b). Reference ideal fringe direction map is depicted in Fig. 6(c).
To sum up the region matching algorithm transforms the preliminary orientation map obtained using PCA into a direction map
Proposed fringe direction estimator employs the Principal Component Analysis method for eigenvectors distribution determination. The PCA method was applied for calculating directions orthogonal and parallel to fringes firstly in . This information was utilized recently to enable efficient ESPI fringe skeletoning , therefore in the following part of the paper we will refer to this algorithm as skeletoning. Our algorithm shows very significant advancement over the one proposed in  - a novel and profitable conditioned region matching algorithm. Using proposed scheme the fringe direction map is calculated for every pixel in the image domain without path-following procedure. No masked out regions encompassing zero-phase-gradient spots are left undetermined (in contrary to skeletoning  where significant number of pixels is undetermined). Additionally there is no risk of error propagation (due to wrong path-following). Moreover the number of algorithm parameters is reduced from three to one simplifying and automating the method. There is no need to set two threshold values due to fully automatic mask designing process. The only parameter is a radius of disk-shaped window for the gradient covariance analysis. Default value is 10 pixels and it enables efficient analysis of a wide variety of fringe patterns in terms of both significant fringe local spacing and shape fluctuations and their defects (noise, contrast and background variations). Algorithm proposed in  provides erroneous results in cases of complicated fringes, like the one in Fig. 1(a). Fringe direction map obtained by scheme reported in  after first unwrapping step is presented in Fig. 7(a) . Higher threshold value determining the result of relative eigenvalue difference map segmentation (Fig. 1(d)) is selected enabling all regions encompassing the zero-phase-gradient spots to be masked out. If the unwrapping process path of the algorithm proposed in  goes through areas marked with blue ellipses the direction angle error appears and starts to propagate. It is caused by very fast fringe direction changes within those areas (see Fig. 1(a)). Problematic areas can be masked out providing sufficiently high threshold value, Fig. 7(b), which results in very unsatisfactory small area of unwrapped direction, Fig. 7(c).
We have implemented second algorithm for comparison purpose - popular combined gradient and plane-fit method proposed in . It generates a fringe orientation map, Fig. 7(d), which has to undergo an unwrapping procedure. We propose to use regular phase unwrapping algorithm  for this purpose. As it is sensitive to 2π jumps the local fringe orientation map should be multiplied by 2 prior to the unwrapping and divided by 2 right after unwrapping. Final direction map, Fig. 7(e), is obtained by 2π wrapping, e.g.,
Proposed direction estimator can be treated as a stand-alone robust algorithm. In our HHT-PCA approach it is aided by modified automatic selective reconstruction fringe filtering and provides the fringe direction map for efficient phase demodulation using the Hilbert spiral transform.
2.2 Modified automatic selective reconstruction algorithm for adaptive fringe pattern filtering
Fringe pattern defects should be removed prior to phase demodulation resulting in augmented accuracy. Fringe pattern description can be simplified to
Proposed by Larkin in [32,33 ] the Hilbert spiral transform is able to calculate amplitude modulation distribution b and phase distribution φ once only the oscillatory term is provided
For fringe pre-filtering - denoising, background removal and contrast enhancement we propose to use a modified automatic selective reconstruction algorithm. The original ASR was recently proposed in  and compared with other popular techniques. Simulations and experimental fringe data analysis corroborated the robustness to noise (additive and multiplicative), significantly varying fringe background, contrast, spatial frequency and orientation (shape). Image is decomposed using very fast EFEMD algorithm into a set of m bidimensional intrinsic mode functions (BIMF) and residue (res), i.e.,
In ASR the filtered image is reconstructed stitching sharply extracted fringe regions from each informative empirical mode - BIMFs from k + 1 to l. k is set automatically to 1 when noise-to-signal ratio of first BIMF is larger than 1; l is set automatically to m. In other words residue is treated as background and first mode (BIMF1) as a noise to be removed. Other modes are considered as informative in some regions. Good regions in each mode are separated using amplitude modulation (b - calculated for each mode by the Hilbert spiral transform using Eq. (11)) value as an indicator. High amplitude corresponds to sharply extracted fringes of interest. For each pixel the mode with highest modulation (BIMFHM) value in this pixel is automatically determined and the corresponding intensity is transferred to reconstructed, filtered fringe pattern
Though ASR performs very efficiently its errors produced by the stitching scheme can be observed as visible stitched region borders . In this contribution we report simple yet effective ASR modification (ASRm) consisting in continuous grayscale weighted analysis. In this approach each pixel from each informative mode contributes to the reconstructed filtered fringe pattern according to its amplitude modulation distribution. Mode amplitudes are calculated using the Hilbert spiral transform and smoothed . Automatic selective reconstruction result is obtained adding informative modes superimposed multiplicatively with their smoothed amplitude modulation distribution (considered as the local fringe quality indicator):28]. Moreover we introduced an option to discard first two modes when noise encountered in fringe pattern is heavy (start with i = 3).
Applying the modified ASR algorithm a considerable fringe pattern quality enhancement is reported. Smoother fringes are obtained in shortened time comparing with the original ASR method. Calculations results corroborating those statements will be presented in the next section.
2.3 Proposed processing path
Proposed all-in-one single-frame HHT-PCA algorithm consists in four steps and provides remarkably accurate phase demodulation in a robust, fast and automatic manner.
- 1. Fringe pattern adaptive pre-filtering by the modified ASR algorithm.
- 2. Filtered fringe direction estimation using PCA based conditioned region matching scheme.
- 3. Filtered fringe efficient phase demodulation using the Hilbert spiral transform aided by the fringe direction map.
- 4. Post-processing by adaptive residual noise reduction.
Steps number 1,2 and 4 are considered as main novel contributions of this paper. In 1 the novelty lies in a parameter-free modified scheme for automatic selective reconstruction of fringe pattern. Errors encountered in the mode stitching process are avoided taking into account for reconstruction every pixel of each informative mode with appropriate weight in the form of smoothed amplitude modulation value. This approach enhances qualitatively and quantitatively the fringe pattern pre-filtering, discards two parameters needed in previous version and disambiguates the final outcome (which previously depended on the parameter choice). In 2 the novelty lies in efficient region matching algorithm which is fast and provides accurate fringe direction value in every pixel - previously reported techniques employing PCA  had to mask out the problematic regions (like phase function extrema or zero phase gradient). In 4 the novelty lies in adaptive removal of residual noise using the EFEMD decomposition of unwrapped phase function and its first mode removal.
In our studies we use advanced and elaborated techniques for fringe pattern pre-filtering which are based on the original 1D Empirical Mode Decomposition (EMD) and apply 2D Hilbert spiral transform which is based on the original 1D Hilbert transform (HT) - this is the reason why we name our multi-step algorithmic solutions as Hilbert-Huang transform (HHT) honoring the original 1D HHT (EMD + HT) scheme .
3. Algorithm evaluation using synthetic fringe data
First step of the proposed HHT-PCA algorithm is the fringe pattern quality enhancement. Proposed novel modified automatic selective reconstruction scheme (ASRm) will be compared with previously reported classical ASR technique and conventional Gaussian filtering followed by the Hilbert transform contrast normalization [25–33 ]. Second step is the fringe direction estimation - proposed PCA algorithm will be compared with the skeletoning technique  and the combined gradient and plane-fit method . Phase demodulation will be conducted using Hilbert spiral transform aided by the direction map. In the proposed HHT-PCA method additional phase filtering in post-processing is implemented. In this way reported all-in-one HHT-PCA algorithm will be compared with two schemes: (A1) Gaussian filtering - fringe direction estimation using the combined method - Hilbert spiral transform phase retrieval and (A2) ASR filtering - fringe direction estimation using skeletoning - HST demodulation. Recently the combination of ASR fringe enhancement with gradient and plane-fit orientation analysis for fringe analysis using Hilbert and Teager transforms was reported in [71, 72 ].
It is important to note that ASR and ASRm procedures are automatic whereas Gaussian filtering needs full supervision. Input fringe pattern, Fig. 1(a), has extremely bad quality resulting in very low quality index QI - a neat quantitative measure of image resemblance compared with the referenced ideal one (proposed by Wang and Bovik in ). Perfect similarity is described by QI = 1. Figure 1(a) has QI = −0.01 which quantifies the extremely cumbersome nature of the simulated pattern generously spoiled by numerous defects.
Quality enhancement results using Gaussian, ASR and ASRm approaches are depicted in Figs. 8(a)-8(c) . Improvement made by the proposed ASRm algorithm over existing ones is readily noticeable, especially in terms of better denoising in regions with lower fringe contrast and no visible border lines. Augmented fringe pattern quality naturally transfers into more accurate local fringe direction estimation using proposed PCA based approach, Fig. 8(f). Previously reported fringe direction estimation schemes like skeletoning  and combined method  provide rather erroneous outcomes regardless the fringe pre-filtering, Fig. 8(d) and 8(e) (due to complicated phase function or intrinsic limitations). Proposed PCA techniques needs around 5 seconds for direction map determination (in current simple, not-optimized and not-parallelized implementation), whereas skeletoning algorithm operates in 19s and the combined method finishes in very fast 0.4s. Calculation time of the latter method greatly augments taking into account the manual parameter setting routine.
It is worth to note the somehow nested enhancement strategy provided by the HHT-PCA algorithm - enhanced fringes result in more accurate fringe direction estimation which ensures higher phase demodulation precision (also provided by the fringe quality improvement). Adaptive filtering in post-processing additionally improves phase quality. Each step contributes to make the proposed single-frame algorithm very robust and attractive.
Phase demodulation results are shown in Fig. 9 . The need for fringe direction map calculation is briefly emphasized presenting wrapped phase map obtained after novel efficient ASRm filtering and Hilbert spiral transform demodulation without any fringe angle information, Fig. 9(a). Ideal simulated wrapped phase is depicted in Fig. 9(b). Proposed technique, Fig. 9(c), show superiority over A1 and A2 - Figs. 9(d) and 9(e), respectively. Unwrapped phase with linear (carrier) term removed for A1 is shown in Fig. 9(f). We decided to not introduce the unwrapped phase obtained by A2 as it is full of erroneous phase jumps produced by incorrect fringe direction map estimation. The final outcome of the proposed processing path, Fig. 9(g), compares favorably with other results in terms of both visual qualitative inspection and QI/RMS computed with reference to ideal phase, Fig. 9(h). Additional adaptive filtering in post-processing in terms of phase distribution EFEMD decomposition and noisy modes removal improves the phase demodulation accuracy (in Fig. 9(g) unwrapped phase after filtering is depicted; in Fig. 9(c) we have presented rewrapped phase distribution after post-filtering). Error map of the best algorithm - the HHT-PCA - is shown in Fig. 9(i).
Similar calculations were performed for the second synthetic fringe pattern, Fig. 10(a) . This time fringes are less complicated to enable successful fringe direction map estimation using previously reported techniques (skeletoning and combined method). Proposed HHT-PCA hybrid algorithm provides better fringe enhancement efficiency and produces more accurate phase distribution than A1 and A2 methods, Fig. 10(b)-10(g), both in terms of quality index and root-mean-square error. Phase error maps computed using simulated reference phase function values show in localized manner the supremacy of the proposed method. Quantitative evaluation using two simulated fringe patterns is summed up in Table 1 corroborating the qualitative visual inspection.
4. Experimental verification
Conclusions derived in the simulation sections are corroborated using experimental fringe patterns. In previous section we have computer-generated ideal fringes, phase and local direction map and treated them as a ground truth for algorithm outcome evaluation (in terms of QI and RMS). In experiments there is no such concept as a ground truth. Single-frame algorithm results will be compared with the reference phase obtained using advanced iterative algorithm  with 9 or 5 phase shifted frames (with pre-filtering).
Fringe analysis results for closed fringes are shown in Figs. 11 and 12 . First fringe pattern has single zero-phase-gradient point (“orientation vortex”) and high spatial frequency fringes which makes it easier to process than the second one (two orientation vortices and low spatial frequency fringes). Quantitative evaluation is presented in Table 2 . Robustness and accuracy of the proposed approach is corroborated. It outperforms previously reported techniques constituting both A1 and A2 schemes.
5. Discussion and conclusions
We have tested using numerical and real experiments three single-frame all-in-one fringe pattern phase demodulation schemes. Two hybrid techniques (A1 and A2) were built innovatively merging already existing algorithmic solutions (A1 - Gaussian filtering, Hilbert normalization, combined gradient and plane-fit for orientation estimation with additional unwrapping to direction map, Hilbert spiral transform phase retrieval; A2 - automatic selective reconstruction for adaptive fringe filtering and normalization, skeletoning for fringe direction estimation and Hilbert spiral transform phase retrieval) while the third one is a completely novel contribution. The proposed HHT-PCA (Hilbert-Huang transform - Principal Component Analysis) employs novel modified automatic selective reconstruction (ASRm) for efficient fringe enhancement. The improvement over classical ASR is made in terms of accuracy boost (region borders influence reduction) and calculation acceleration (due to automatic background identification procedure and possible first modes neglecting basing on a priori knowledge). Hilbert spiral transform is used for phase demodulation; additional adaptive EFEMD filtering of phase distribution reduces residual noise transferred from fringe domain to phase domain. Therefore the general concept of Hilbert-Huang transform (HHT) is responsible for fringe filtering and demodulation. The latter would not be possible without robust and accurate fringe direction map estimation. For this reason a Principal Component Analysis based novel conditioned region matching technique was implemented. The improvement over existing solutions is made in terms of automatic operation with no orientation unwrapping needed (comparing with manually supervised combined gradient and plane-fit method followed by unwrapping) and enhanced accuracy full-field direction estimation (comparing with the path-following, heavily masked skeletoning technique; very susceptible to multiple phase local extrema - “orientation vortices”).
Proposed HHT-PCA technique was proven uniquely robust and extremely effective for a wide range of complicated fringe patterns with noise, rapid contrast changes, uneven background influence, significant local fringe orientation and period variations. Comprehensive qualitative and quantitative studies have shown impressive accuracy, up to RMS = 0.03 radians, of the HHT-PCA - it compares very favorably with other existing cutting-edge single-frame techniques. It is very important to note that it also outperforms popular two-frame methods like Gram-Schmidt orthonormalization with HHT pre-filtering  or Optical Flow , with RMS up to 0.2-0.6 radians. Excellent features of the proposed HHT-PCA method follow from its global-local properties (global Hilbert transform, pixel-by-pixel ASRm scheme and region-by-region PCA direction estimator) and nested improvement strategy (each step of the processing path improves the outcome of the other steps).
We advocate the use of the HHT-PCA all-in-one single-frame phase retrieval algorithm because of its superb and unparalleled robustness, versatility and accuracy corroborated by successful analysis of synthetic and experimental complicated fringe patterns. It can be successfully applied to, e.g., optical element testing, experimental mechanics or life sciences  aiding various fringe-based optical measurement techniques.
This work was supported, in part, by National Science Center, Poland, under the projects OPUS 8 (UMO-2014/15/B/ST7/04650) and PRELUDIUM 8 (UMO-2014/15/N/ST7/04881), in part by statutory funds and grant founded by the Dean of Faculty of Mechatronics, Warsaw University of Technology. Work of MT was partially supported by the Polish Ministry of Science and Higher Education and Foundation for Polish Science. Authors thank Z. Sunderland for providing experimental data used in Figs. 11 and 12.
References and links
1. J. Schwider, “Advanced evaluation techniques in interferometry,” in Progress in Optics, E. Wolf ed. (Elsevier, 1990).
2. D. W. Robinson and G. Reid, Interferogram Analysis: Digital Fringe Pattern Measurement (Institute of Physics Publishing, 1993).
3. D. Malacara, ed., Optical Shop Testing (Jon Wiley and Sons, 2007).
4. M. Servin, J. A. Quiroga, and M. Padilla, Fringe Pattern Analysis for Optical Metrology: Theory, Algorithms, and Applications (Wiley, 2014).
5. M. Takeda, H. Ina, and S. Kobayashi, “Fourier-transform method of fringe-pattern analysis for computer-based topography and interferometry,” J. Opt. Soc. Am. 72(1), 156–160 (1982). [CrossRef]
6. T. Kreis, “Digital holographic interference-phase measurement using the Fourier-transform method,” J. Opt. Soc. Am. A 3(6), 847–855 (1986). [CrossRef]
7. J. H. Massig and J. Heppner, “Fringe-pattern analysis with high accuracy by use of the Fourier-transform method: theory and experimental tests,” Appl. Opt. 40(13), 2081–2088 (2001). [CrossRef] [PubMed]
9. Q. Kemao, “Two-dimensional windowed Fourier transform for fringe pattern analysis: Principles, applications and implementations,” Opt. Lasers Eng. 45(2), 304–317 (2007). [CrossRef]
10. Q. Kemao, Windowed Fringe Pattern Analysis (SPIE, 2013).
12. K. Pokorski and K. Patorski, “Visualization of additive-type moiré and time-average fringe patterns using the continuous wavelet transform,” Appl. Opt. 49(19), 3640–3651 (2010). [CrossRef] [PubMed]
14. 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]
15. M. Servin, J. L. Marroquin, and F. J. Cuevas, “Fringe-follower regularized phase tracker for demodulation of closed-fringe interferograms,” J. Opt. Soc. Am. A 18(3), 689–695 (2001). [CrossRef]
17. M. Kujawinska and J. Wojciak, “Spatial-carrier phase-shifting technique of fringe pattern analysis,” Proc. SPIE 1508, 61–67 (1991). [CrossRef]
18. M. Pirga and M. Kujawinska, “Two directional spatial-carrier phase-shifting method for analysis of crossed and closed fringe patterns,” Opt. Eng. 34(8), 2459–2466 (1995). [CrossRef]
19. K. Creath and J. Schmit, “N-point spatial phase-measurement techniques for non-destructive testing,” Opt. Lasers Eng. 24(5–6), 365–379 (1996). [CrossRef]
20. A. Styk and K. Patorski, “Analysis of systematic errors in spatial carrier phase shifting applied to interferogram intensity modulation determination,” Appl. Opt. 46(21), 4613–4624 (2007).
22. Y. Du, G. Feng, H. Li, J. Vargas, and S. Zhou, “Spatial carrier phase-shifting algorithm based on principal component analysis method,” Opt. Express 20(15), 16471–16479 (2012). [CrossRef]
24. X. Zhu, C. Tang, B. Li, C. Sun, and L. Wang, “Phase retrieval for single frame projection fringe pattern with variational image decomposition,” Opt. Lasers Eng. 59(8), 25–33 (2014). [CrossRef]
26. M. Trusiak and K. Patorski, “Space domain interpretation of incoherent moiré superimpositions using FABEMD,” Proc. SPIE 8697, 869704 (2012). [CrossRef]
27. 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]
29. 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]
30. K. Patorski, M. Trusiak, and M. Wielgus, “Fast adaptive processing of low quality fringe patterns by automated selective reconstruction and enhanced fast empirical mode decomposition,” Fringe 2013, 185–190 (2014).
31. M. Trusiak, K. Patorski, and M. Wielgus, “Hilbert-Huang processing and analysis of complex fringe patterns,” Proc. SPIE 9203, 92030K (2014). [CrossRef]
32. K. G. Larkin, D. J. Bone, and M. A. Oldfield, “Natural demodulation of two-dimensional fringe patterns. I. General background of the spiral phase quadrature transform,” J. Opt. Soc. Am. A 18(8), 1862–1870 (2001). [CrossRef] [PubMed]
33. K. G. Larkin, D. J. Bone, and M. A. Oldfield, “Natural demodulation of two-dimensional fringe patterns. II. Stationary phase analysis of the spiral phase quadrature transform,” J. Opt. Soc. Am. A 18(8), 1871–1881 (2001). [CrossRef] [PubMed]
35. B. Jahne, Practical Handbook on Image Processing for Scientific Applications (CRC, 1997).
38. Q. Yu and K. Andresen, “Fringe-orientation maps and fringe skeleton extraction by the two-dimensional derivative-sign binary-fringe method,” Appl. Opt. 33(29), 6873–6878 (1994). [CrossRef] [PubMed]
40. Q. Yu, X. Liu, and X. Sun, “Generalized spin filtering and an improved derivative-sign binary image method for the extraction of fringe skeletons,” Appl. Opt. 37(20), 4504–4509 (1998). [CrossRef] [PubMed]
41. J. L. Marroquin, R. Rodriguez-Vera, and M. Servin, “Local phase from local orientation by solution of a sequence of linear systems,” J. Opt. Soc. Am. A 15(6), 1536–1544 (1998). [CrossRef]
45. Q. Yu, X. Yang, S. Fu, and X. Sun, “Two improved algorithms with which to obtain contoured windows for fringe patterns generated by electronic speckle-pattern interferometry,” Appl. Opt. 44(33), 7050–7054 (2005). [CrossRef] [PubMed]
46. X. Yang, Q. Yu, and S. Fu, “An algorithm for estimating both fringe orientation and fringe density,” Opt. Commun. 274(2), 286–292 (2007). [CrossRef]
47. X. Yang, Q. Yu, and S. Fu, “A combined method for obtaining fringe orientations of ESPI,” Opt. Commun. 273(1), 60–66 (2007). [CrossRef]
48. X. Yang, Q. Yu, and S. Fu, “Determination of skeleton and sign map for phase obtaining from a single ESPI image,” Opt. Commun. 282(12), 2301–2306 (2009). [CrossRef]
49. F. Zhang, W. Liu, J. Wang, Y. Zhu, and L. Xia, “Anisotropic partial differential equation noise-reduction algorithm based on fringe feature for ESPI,” Opt. Commun. 282(12), 2318–2326 (2009). [CrossRef]
50. A. M. Siddiolo and L. D’Acquisto, “A direction/orientation- based method for shape measurement by shadow Moire,” IEEE Trans. Instrum. Meas. 57(4), 843–849 (2008). [CrossRef]
51. C. Tang, L. Han, H. Ren, D. Zhou, Y. Chang, X. Wang, and X. Cui, “Second-order oriented partial-differential equations for denoising in electronic-speckle-pattern interferometry fringes,” Opt. Lett. 33(19), 2179–2181 (2008). [CrossRef] [PubMed]
54. C. Tang, Z. Wang, L. Wang, J. Wu, T. Gao, and S. Yan, “Estimation of fringe orientation for optical fringe patterns with poor quality based on Fourier transform,” Appl. Opt. 49(4), 554–561 (2010). [CrossRef] [PubMed]
55. H. Wang, Q. Kemao, R. Liang, H. Wang, M. Zhao, and X. He, “Oriented boundary padding for iterative and oriented fringe pattern denoising techniques,” Sig. Proc. 102(9), 112–121 (2014). [CrossRef]
56. D. Zhang, M. Ma, and D. D. Arola, “Fringe skeletonizing using an improved derivative sign binary method,” Opt. Lasers Eng. 37(1), 51–62 (2002). [CrossRef]
57. C. Tang, W. Lu, Y. Cai, L. Han, and G. Wang, “Nearly preprocessing-free method for skeletonization of gray-scale electronic speckle pattern interferometry fringe patterns via partial differential equations,” Opt. Lett. 33(2), 183–185 (2008). [CrossRef] [PubMed]
61. M. Servin, J. A. Quiroga, and J. L. Marroquin, “General n-dimensional quadrature transform and its application to interferogram demodulation,” J. Opt. Soc. Am. A 20(5), 925–934 (2003). [CrossRef] [PubMed]
62. D. Crespo, J. A. Quiroga, and J. A. Gomez-Pedrero, “Fast algorithm for estimation of the orientation term of a general quadrature transform with application to demodulation of an n-dimensional fringe pattern,” Appl. Opt. 43(33), 6139–6146 (2004). [CrossRef] [PubMed]
63. J. A. Quiroga, M. Servin, J. L. Marroquin, and D. Crespo, “Estimation of the orientation term of the general quadrature transform from a single n-dimensional fringe pattern,” J. Opt. Soc. Am. A 22(3), 439–444 (2005). [CrossRef] [PubMed]
64. J. Villa, I. De la Rosa, G. Miramontes, and J. A. Quiroga, “Phase recovery from a single fringe pattern using an orientational vector-field-regularized estimator,” J. Opt. Soc. Am. A 22(12), 2766–2773 (2005). [CrossRef] [PubMed]
66. J. C. Estrada, M. Servín, J. A. Quiroga, and J. L. Marroquín, “Path independent demodulation method for single image interferograms with closed fringes within the function space C2.,” Opt. Express 14(21), 9687–9698 (2006). [CrossRef] [PubMed]
71. B. Deepan, C. Quan, and C. J. Tay, “Phase retrieval from a single fringe pattern using Teager-Hilbert- Huang transform,” Proc. SPIE 9302, 93020Q (2015). [CrossRef]
72. B. Deepan, C. Quan, and C. J. Tay, “Determination of phase derivatives from a single fringe pattern using Teager-Hilbert-Huang transform,” Opt. Commun. 359, 162–170 (2016). [CrossRef]
74. H. Wang and Q. Kemao, “Quality-guided orientation unwrapping for fringe direction estimation,” Appl. Opt. 51(4), 413–421 (2012). [PubMed]
76. J. A. Quiroga, M. Servin, and F. Cuevas, “Modulo 2π fringe orientation angle estimation by phase unwrapping with a regularized phase tracking algorithm,” J. Opt. Soc. Am. A 19(8), 1524–1531 (2002). [CrossRef] [PubMed]
77. X. Zhou, J. P. Baird, and J. F. Arnold, “Fringe-orientation estimation by use of a Gaussian gradient filter and neighboring-direction averaging,” Appl. Opt. 38(5), 795–804 (1999). [CrossRef] [PubMed]
78. J. Vargas, J. A. Quiroga, C. O. S. Sorzano, J. C. Estrada, and J. M. Carazo, “Two-step interferometry by a regularized optical flow algorithm,” Opt. Lett. 36(17), 3485–3487 (2011). [CrossRef] [PubMed]
79. M. A. Herráez, D. R. Burton, M. J. Lalor, and M. A. Gdeisat, “Fast two-dimensional phase-unwrapping algorithm based on sorting by reliability following a noncontinuous path,” Appl. Opt. 41(35), 7437–7444 (2002). [CrossRef] [PubMed]
80. N. E. Huang, Z. Sheng, S. R. Long, M. C. Wu, W. H. Shih, Q. Zeng, N. C. Yen, C. C. Tung, and H. H. Liu, “The empirical mode decomposition and the Hilbert spectrum for non-linear and non-stationary time series analysis,” Proc. R. Soc. Lond. A 454(1971), 903–995 (1998). [CrossRef]
81. Z. Wang and A. C. Bovik, “A universal image quality index,” IEEE Signal Process. Lett. 9(3), 81–84 (2002). [CrossRef]
83. M. Trusiak and K. Patorski, “Two-shot fringe pattern phase-amplitude demodulation using Gram-Schmidt orthonormalization with Hilbert-Huang pre-filtering,” Opt. Express 23(4), 4672–4690 (2015). [CrossRef] [PubMed]