Abstract
Due to the physical nature of the interference phenomenon, extracting the phase of an interferogram is a known sinusoidal modulation problem. In order to solve this problem, a new hybrid mathematical optimization model for phase extraction is established. The combination of frequency guide sequential demodulation and harmony search optimization algorithms is used for demodulating closed fringes patterns in order to find the phase of interferogram applications. The proposed algorithm is tested in four sets of different synthetic interferograms, finding a range of average relative error in phase reconstructions of 0.14–0.39 rad. For reference, experimental results are compared with the genetic algorithm optimization technique, obtaining a reduction in the error up to 0.1448 rad. Finally, the proposed algorithm is compared with a very known demodulation algorithm, using a real interferogram, obtaining a relative error of 1.561 rad. Results are shown in patterns with complex fringes distribution.
1. Introduction
Interferometry includes a series of techniques used in the measurement of aberrations, deformations, flatness, and perturbations. It may also be applied to measure variables as temperature gradients, strain analysis, depth measurement, and so forth. Widely known techniques to demodulate interferograms are phase shifting algorithms for several images and the Fourier method for a single linear carrier fringe pattern [1]. However, these techniques are difficult to implement when the object under study changes fast and continuously, and the dynamic range of the phase does not allow the use of a large linear carrier without infringing the sampling theorem. In those situations a single interferogram with closed and possibly complex fringe distribution must be analyzed to recover the phase information related to the physical phenomena being measured. This is known as a difficult task since there are many solutions that are compatible with the measured data but lacks physical meaning. The accuracy of measurements carried out from a single fringe pattern that includes closed fringes is thus intensely dependent on the phase distribution of the recorded interferogram being estimated. Recently, many phase recovering methods have been developed as combination of genetic algorithms and parametric methods [2–4], soft computing techniques applied to Zernike polynomials [5], combination of genetic algorithms and frequency guided sequential demodulation [6], particle swarm optimization [7], unwrapping of phase maps with sign changes [8], two-dimensional regularized phase-tracking technique [9], and so forth. In general, there are not particular processes that succeed in obtaining the phase for any given interferogram, but all of them are limited to specific features of the fringe pattern.
The phase demodulation trouble has been formulated as an optimization challenge, where soft computing procedures may be used to find the phase solution that best matches the nonlinear equation represented by fringe patterns. Few years ago, genetic algorithms have been tested [2]; the authors developed a parametric method for fringe pattern demodulation using a genetic algorithm (GA). A parametric estimation of the coefficients of a 15th degree Zernike polynomial is used in order to approximate the phase; a population of chromosomes is programmed within the coefficients to calculate the phase. A cost function is then employed considering the number of the observed fringes and the fringes that result from the recovered phase match, the phase softness, and the prior knowledge of the object. Normally, the final solution of the GA is based on a cost function, which is stated as the comparison between the better individual in the population and the target (real fringes); a population evolution process is allowed until a cost function average threshold is achieved. The authors reported a root mean square (rms) error of 0.12 radians. This method was applied to noisy fringe patterns and to a single closed fringe image. Additional improvements and variations of this work were subsequently presented by the same research team [3, 4].
Another soft computing technique used for phase reconstruction is particle swarm optimization (PSO). This algorithm was introduced by Kennedy and Eberhart in 1995 [10], as an evolving optimization technique. In 2012 Jiménez et al. [7] used PSO for phase recovery; they compared a GA and a PSO for phase recovery on several fringe patterns, obtaining errors of 0.4281 and 0.313 rad., respectively, showing an improvement in accuracy of PSO over GA; processing time improvements were announced, but no results were shown.
As mentioned before, the demodulation of a single interferogram often involves a combination of methods (GA + Zernike, PSO + Zernike, Neuronal networks + others, etc.). In 2009 Wang and Kemao reported a new hybrid method; they used frequency guided sequential demodulation (FSD) as interferogram demodulator, combined with Levenberg Marquardt (LM) optimization [11], method implemented by their quickness and efficiency in fringes demodulation.
In this work a FSD with harmony search optimization (HSO) is investigated in order to test the performance in a single interferogram with closed fringes. The main motivation is the advantages of the HSO technique over other soft computing techniques already reported. The HSO technique was inspired in the observation of musical composition to search a perfect harmony and was introduced by Geem et al. in 2001 [12] and has found its way in several applications as diverse as engineering, math, industrial process, biology, and so forth [13–19]. An excellent recent review and categorization of the applications of HSO was conducted by Manjarres et al. in 2013 [20]. Some advantages of this method are that it uses simple algebraic equations and real values, while the derivative information is unnecessary unlike GA and other optimization techniques.
In the following section, the physical theory of the interferograms is presented as well as the concepts of the HSO and FSD algorithms. In the next section, the image-processing techniques and the experimental setup used to implement the soft computing proposed method are described. Finally, in the last two sections the results and the conclusions are presented, respectively.
2. Theory
Metrology has techniques such as fringe projection profilometry and optical interferometry to measure physical quantities in many areas of engineering and science, but the importance of these methods lies in the fact that they are noninvasive procedures [21]. Recently, advances in computational techniques have the potential to extend the measuring capabilities of optical metrology applications. In the present section the optical metrology basis, the harmony search optimization model, and the frequency guided sequential demodulation process are shown.
2.1. Optical Interferometry
Interferometry studies the engagement of two or more light waves, where one of them has suffered a modification by one characteristic of an object being tested [22]. Demodulation of the phase is the most important task in interferometry measurements; the phase is related to a physical quantity to be measured.
The optical arrangement, shown in Figure 1, is a Twyman-Green interferometer setup (a Michaelson interferometer modification). The interference is produced by the difference of optical path between the two arms of the interferometer. The interferogram is reordered by a photo detector array (e.g., a charged coupled display camera) and then digitized for show on a monitor or stored for further processing with computational algorithms like unwrapped phase, digital filtering, demodulation phase, and so forth [23]. The optical components Lc, Lf, and Bs are a positive collimating lens, a positive focusing lens, and a beam splinter, respectively. The fringe pattern intensity is modeled bywhere , and are the interference fringe pattern intensity, the background illumination, the modulation amplitude, the phase term, and the spatial coordinates of the surface under test, respectively.
2.2. Harmony Search Optimization
The harmony search is a metaheuristic optimization algorithm; it was created by Lee and Geem in 2005 [24]; it was inspired in music harmony and is a powerfully soft computing tool to solve many optimization problems. The HSO algorithm is structured in five steps [25], as shown in Figure 2. Those steps are the following: initialize the problem and algorithm parameters; initialize the Harmony Memory (HM) array; improvise a new harmony from HM; update the HM; and check the stopping criterion. These steps are described as follows.
Step 1. The optimization problem is specified and initialized its parameter values as follows:where , and are the objective function, the set of each decision variable , and the set of the possible range of values for each decision variable, respectively. Also, the search range of is , with and as the lower and upper bounds for each decision variable, respectively. The parameters to initialize the algorithm are the number of decision variables (NVAR); the number of improvisations or stopping criterion (SC); the number of solution vectors in the Harmony Memory or the Harmony Memory Size (HMS); Harmony Memory Considering Rate (HMCR); Pitch Adjusting Rate (PAR) for each generation; and the arbitrary distance bandwidth for each generation (). In the following steps the role of some of these parameters is explained.
Step 2. The Harmony Memory (HM) array is complete with a random solution vectors as the HMS, as can be seen in the following:where HM acts as a dynamical memory address, where the sets of decision variables (all the solution vectors) are saved. This HM is similar to the genetic array in the genetic algorithm (GA) [26].
Step 3 (improvise a new harmony with parameters adjustment). For this improvisation, a new harmony vector is generated taking into account three basic considerations [27]: memory consideration, pitch adjustment, and random selection.
In the memory consideration, the value of the first decision variable for the new vector is chosen from any of the values in the specified HM. The values of other decision variables are chosen in the same manner. The HMCR has a range of values from 0 and 1, while is the possibility of creating a random HM member [27].
In the second consideration, a pitch adjusting is performed only after a value is chosen from the HM. The PAR considers the possibility of change of some elements of HM.
In the last consideration, a new vector is generated with a pitch adjustment threshold of 0.5, and the definition of this vector is described bywhere has an uniform distribution between −1 and 1. The function is fed with a random number between 0 and 1 to adjust the PAR. The sign of the addition in (4) depends of the value obtained with function; if the value is less than 0.5 the sign is positive, and negative in other case.
Other parameters used in this step are NV and GEN, a new vector to improvise HM and for the depuration of the actual generation, respectively.
Step 4 (update HM). If the new harmony vector is better than the worst harmony vector in HM, then this latter is replaced by the new harmony vector. The HM is then sorted by the objective function value.
Step 5 (check the stopping criterion). If the SC is satisfied, computation is terminated. Otherwise, Steps 3 and 4 are repeated.
The HMCR and PAR parameters introduced in Step 3 help the algorithm to find locally and completed improved solutions, respectively. PAR and in HS algorithm are important parameters in fine-tuning of the optimized solution vectors and can be potentially useful in adjusting convergence rate of the algorithm to an optimal solution. The HSO pseudocode is shown in Algorithm 1.
| ||||||||||||||||||||||||||||||
The traditional HS algorithm uses fixed value for both PAR and . In the HS method PAR and values adjusted in initialization step (Step 1) and cannot be changed during new generations. Large values with small PAR values can produce poor performance of the algorithm and increased iterations are needed to find an optimum solution. Although in early generations must take a bigger value to enforce the algorithm to increase the diversity of solution vectors, small values in final generations increase the fine-tuning of solution vectors. Also, small values with large PAR values usually cause the improvement of best solutions in final generations in which the algorithm converges to the optimal solution vector [28].
2.3. Frequency Guided Sequential Demodulation
The frequency guided sequential demodulation (FSD) method was created by Kemao and Soon in 2007 [29], and it is used for recovery from the phase of closed fringe patterns. The algorithm of this method is explained in the following six sentences.
The first sentence is fringe denoising. The noise of (1) will be removed (low-pass filtering).
The second sentence is normalization. The parameters of resulting closed interferogram are normalized [30]; this is the background and the amplitude ; in effect, (1) is modified in order to find the interference fringe pattern normalized intensity
The third sentence is extracting the phase. The direct phase is extracted by the following equation:which is used to calculate the iteration frequencies or intermediate frequencies and ; those frequencies are temporal, whilst they refine optimization in order to find the final frequencies and ; this is, ; it is feasible to use exhaustive search algorithm that guarantees a global minimum [31].
The fourth sentence is about the extraction of local frequency. The true phase is to be considered locally linear, aswhere and the operator “” are the coordinates of the pixel under study, neighbors pixel coordinates, the internal frequencies while the optimization process is being performed, and the dot product, respectively. In other words is the intermediate parameter vector. Then, a virtual fringe pattern can be generated as
The fitness or energy cost function to be used by the HSO algorithm is equivalent to the squared difference between the virtual and the real fringe pattern; that is,where is a subimage window function briefly studied. The local frequency is estimated by minimizing the fitness function:
The fifth sentence is the extraction of estimated frequencies. Finally, the two coefficients found for the pixel under study or the estimated frequencies by HSO are
The sixth sentence is eliminated the ambiguity error. In order to determine the frequency guided sign, all subimages were processed to correct the problem of ambiguity sign. Palpably and have sign ambiguity and the sign function can be determined by forcing the local frequencies to be continuous; this is,and the ambiguity is indicated by the subindex “.” The sign determination is then continued with a pixel that adjoins pixel until all of the pixels have been processed. Once the sign field is determined, both and can be determined [31] by
The Levenberg-Marquardt optimization is another algorithm that has been applied for this kind of recovery phase; this algorithm uses the advantages of Gauss-Newton and gradient-descent methods based on adaptive rules [11].
3. Fringe Pattern Demodulation with FSD-HSO Algorithm
As explained above, (5) expresses the normalized interferogram and (6) represents the direct phase extracted. The process to demodulate a closed fringe pattern applying the FSD-HSO algorithm consists in splitting the complete closed fringe pattern into subimages of 5 × 5 pixels (see Figure 3). Additionally, a pixel is demodulated with the nearest eight neighbors, and the HSO algorithm is implemented after selecting the pixel in order to recover the phase for that element.
Further, the technique FSD is solved applying the HSO algorithm and taking into consideration (7), where the true phase is locally linear. Equation (11) expresses the two calculated frequencies by the HSO algorithm. Equation (8) describes the interferogram obtained with the estimated phase. Equation (9) depicts the fitness function for the HSO; it is the square of the subtraction of the real and the calculated fringe patterns; hence the proper frequencies cause an energy cost function minimization. The above process is applied in all pixels for the subimage. When the process is finished for all the subimages, (13) corrects the ambiguity error. Figure 4 shows the diagram of the demodulation process for each pixel applying HSO and FSD.
The pseudocode algorithms necessary to implement the proposed FSD-HSO procedure are shown in Algorithms 1, 2, 3, 4. First algorithm is the HSO method; second algorithm is the sign determination, in order to find the sign ambiguities of the estimate frequencies; third algorithm implements the procedure to clean the interferograms; fourth algorithm is the local frequencies calculation.
| ||||||||||||||
| ||||||||||||||
| ||||||||||||||||||||||
4. Results of Closed Fringe Pattern Demodulation with FSD-HSO Algorithm
A variety of physical quantities can be measured through the phase; it can be using interferometric techniques. The method to calculate the phase for a closed fringe pattern is explained below. Synthetic interferograms (closed fringe) are demodulated using the FSD-HSO algorithm. The physical variable associated with the obtained phase is smooth and continuous.
The HSO algorithm begins its search with HM and in this case is created by rows of individuals with two elements (NVAR), each with random values that correspond to the beginning of estimated frequencies by (14); this isThe values of each parameter in HSO (NVAR, SC, HMS, PAR, HMCR, and ) are according to the interferogram to be demodulated. In this work, four distinct simulated closed fringe patterns (synthetic interferograms, SIs) and a real interferogram (RI) are exposed. The first example, synthetic-interferogram-one (SI1), corresponds to Newton rings interferogram, that is, classical defocus interference generated in a Twyman-Green or a Fabry-Perot interferometer; the synthetic-interferogram-two (SI2) is a pair of symmetrical positive lobules; the synthetic-interferogram-three (SI3) is a pair of asymmetrical positive lobules; finally, the synthetic-interferogram-four (SI4) is a pair of inverse lobules; RI is a noisy interferogram with defocus aberration. These five examples are a compendium of different difficulties of closed fringe demodulation. The normalization of the SIs and RI were calculated and shown in Figure 5(a). The normalized interferogram was partitioned into subimages of 5 × 5 pixels and also is considering the linear phase to be local. The direct phase is created by inverse cosine of normalized image; see Figure 5(b). Each subimage was processed and estimated the frequencies and , along and directions, respectively; see Figures 5(c) and 5(d). Equation (10) represents the cost function; it is decreased close to the selected pixel and its eight neighbors. A small value in cost function implies good performance in the demodulation. The direct phase was processed by the FSD-HSO algorithm described in flow diagram of Figure 4; as seen in Figure 5(e), the sign ambiguity creates discontinuities in this step of phase reconstruction process. In Figure 4 a block of sign determination was included because of the phase ambiguity problem presented in the function cosine; this is ; details of this method can be found in [31]. Once the subimage is completed, the HSO starts with other subimages; when the last subimage is demodulated, the procedure finishes. At last, the estimated phases are shown in Figure 5(f).
(a)
(b)
(c)
(d)
(e)
(f)
In order to observe minutiae details of the synthetically estimated phases, a series of 3D graphs are presented in Figure 6. The estimated phase of SI1 is observed as a single negative paraboloid, as shown in Figure 6(a); in counterpart, the other three SIs estimated phases have two peaks. The estimated phase of SI2 is a pair of positive symmetrical peaks; see Figure 6(b). However, the respective SI3 estimated phase, is observed as two positive asymmetrical peaks; see Figure 6(c). At last, two inverse peaks are represented in Figure 6(d) as SI4 estimated phase. This last interferogram was the most difficult to demodulate, because the search process of the HSO algorithm was slow due to the increase in iterations for the PAR and adjustment. It is seen that the combination of the HSO and FSD methods can handle interferograms with more complex fringe distribution than circular fringes.
(a)
(b)
(c)
(d)
The average relative error was calculated by the following equation [3]:where , are the original interferogram phase and the interferogram size, respectively.
In addition to the synthetic analysis, the RI estimated phase demodulated with FSD-HSO was compared with a robust unwrapped phase method, implemented by [8], for this occasion that method is taken to calculate the original interferogram phase. In order to compare the two results visually, their 3D phases are shown in Figures 7(a) and 7(b); these two graphs are very similar; to complement the analysis, the variance between the original phase and the calculated phase was obtained, for this RI whose difference was very small (around ~0.5 radians), as shown in Figure 7(c). For SI analysis, the original phase is the phase with which was programmed each synthetic interferogram.
(a)
(b)
(c)
Additionally, the results of HSO algorithm are compared with genetic algorithm (GA) method, both using FSD as demodulation technique. Results of this comparison can be seen in Table 1, where column one contains the five different interferograms used in the present work. The second column of Table 1 contains the peak to valley deviations, where column represents the maximal and minimal values obtained in the respective estimated phases. The third column of the same table presents the rms deviations of the estimated phases. In and rms calculus significant differences cannot be observed between GA and HSO methods, but these data are presented because of their physical significance. The average relative errors were presented in fourth column; it can be observed that HSO algorithm has smaller errors than GA algorithm; the difficulties presented in SI4 and RI (pronounced sign peaks asymmetry and/or very noisy image) increment a little the error in both methods; the biggest error was 1.561 rad.; therefore, the robustness of FSD-HSO is confirmed. The last column presents the time processing, where the FSD-HSO demonstrated as a fast method with symmetrical sign peaks interferograms demodulation (until 1.5 minutes), but the time processing increases as different sign peak lobules appear in the interferogram. Results of Table 1 are supported by the plot of Figure 8, in which the number of iterations for the proposed method is about one half to achieve convergence.
The selected parameters used in this work are shown in Table 2. The entire SIs and RI images had MN size of 250 × 250 pixels and 128 × 128 pixels, respectively. In order to obtain consistent results, all the calculus is made with the same computational system: Intel Celeron B815 1.6 GHz, with 2 GB of RAM, 32-bit operative system, Windows 7, and Matlab R2011a 7.120.635. The FSD-GA at its best operation, used to compare the proposed algorithm, has the following parameters: individuals 1500–2500; bits by gen 12–15; gens 2; generations 1000–1500; limits of each coefficient []–[].
Experimentally it is observed that the proposed algorithm’s tuning parameters depend on the interferogram characteristics. The best resolution in the image demodulation was found experimentally where PAR must approach its upper limit (one) and tends to its lower limit (zero).
5. Conclusions
The use of FSD algorithm to demodulate the phase is a simple technique of soft computing, because just two variables for this work, and , need to be calculated compared to other polynomial approaches that calculate more variables or coefficients.
Instead, the HSO is a relatively new optimization technique that requires a few input parameters. The fitness function used by the HSO is small and the data search is relational, so convergence is fast and has a low percentage of error.
Therefore, the combined algorithm FSD-HSO is an uncomplex algorithm because of its algebraic nature, and, used in the demodulation of closed fringe patterns, its convergence is fast. The error calculated in synthetic interferograms simulated was 0.14–0.39 which is a very good approximation. The method is limited only to demodulation of closed fringe patterns.
The proposed algorithm is relatively fast compared to the FSD-AG algorithm; it converges with less iterations. It only presents difficulties in interferograms with different sign peak lobules, but this disadvantage is present in both methods.
Other methods to demodulation also shed low error rates, such as the AG-Zernike algorithm, but, using these two methods for example, the algorithm becomes very complex, since it is necessary to perform a lot of computational processes such as data mutations, data changes from binary to real, and a lot of variables to be estimated (a minimum of 7 variables).
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
Mr. Rodriguez-Marmolejo would like to acknowledge Centro Universitario de los Lagos and Instituto Tecnológico de Aguascalientes authorities for the support in his Ph.D. studies.