## Abstract

In this paper we present several eight-frame algorithms for their use in phase shifting profilometry and their application for the analysis of semi-fossilized materials. All algorithms are obtained from a set of two-frame algorithms and designed to compensate common errors such as phase shift detuning and bias errors.

© 2013 OSA

## 1. Introduction

The analysis of interferograms can be categorized in two methods: temporal methods (phase shifting) and spatial methods [1]. It is well known that in the temporal phase shifting (TPS) techniques the existing algorithms exhibit considerable measurement inaccuracy unless phase shifts are precisely known. Accordingly, many efforts have been done in recent decades to establish effective error-compensating phase-extraction algorithms [2,3]. Most of the existing phase-shifting algorithms are based on the assumption that the phase-shift applied to each pixel of the intensity frame is a known constant value. However, it may be very difficult to achieve this in practice, mainly because one of the enduring problems with temporal phase-shifting algorithms (PSA) is precisely the phase-shift errors or miscalibration. By using more than three frames, it is possible to design algorithms to compensate the deterministic shift errors (such as nonlinearities of first, second and third order for the piezoelectric), and others like non-sinusoidal fringe profiles [1]. These TPS techniques are applied for the design of non-contact optical 3D-profiling instruments, which are often used for the study of surface modifications of mammalian fossil bones. The main advantage of optical laser scanning of fossil and semi-fossil bones is that the laser’s non-contact nature permits the analysis of a small fragile and poorly preserved surface [4]. Therefore, the laser-stylus 3D-microprofiling has revealed several qualities that make it particularly suitable for the study of fossil and semi-fossil bones, as they cause no alteration to the surface, allowing to study objects under atmospheric conditions. In addition, complex surface features such as flank inclinations, symmetry and micro-striations are accessible for quantitative evaluation, and also features hidden by shallow overall relief are also visually extracted [5]. Considering the importance in different research fields of having a full digital documentation of fossilized bones, for example to compare morphological measurements, it is critical to have a means of obtaining accurate digital three-dimensional images. Moreover, digital three-dimensional imaging becomes crucial for studies that involve the use of deoxyribonucleic acid (ancient DNA) in which samples (semi-fossil in this case) need to be completely pulverized to extract DNA. Well preserved fossil samples, with minimum handling, are essential in order to obtain high quality DNA. This was the starting point for our interest to design phase filters that satisfy the requirements of the optical experimental setup. The technique of fringe projection profilometry (FPP) is used for the recovering of surface topography. Nevertheless, the accuracy of this technique is limited by the presence of different systematic and random error sources, such as higher harmonics in the intensity signal, phase-shifter miscalibration, nonlinear response of the photodetectors, nonsinusoidal periodic waveforms, random intensity noise, speckle decorrelation, and vibration. All algorithms compensate for some kind of error, thus in the experiment we present, we calculated several filters with specific behaviors. As the errors cannot be completely eliminated but only minimized [1], the common technique to minimize systematic errors is the use of a quadrature filter insensitive to such errors [6–8]. There are already several methods to design quadrature filters and a large number of algorithms in the literature [6–13]. However, most of them only give a particular algorithm designed for a very specific phase step, which are not tuned on the ideal phase step. Hence, in this study we introduce several eight-frame filters with phase steps of $\pi /4$ that have not been reported. These algorithms are designed to compensate some systematic errors such as miscalibration and bias variation, which implies more accuracy and better signal recovery for eight step systems. The present work offers a new set of eight-frame algorithms for phase extraction, which are obtained from a set of two-frame filters and are designed to achieve the best signal-to-noise ratio (SNR), capable of minimizing, and even compensate the most of the systematic errors as quadrature and detuning errors. The algorithms are tested and their efficiency corroborated by using computer simulated interferograms and a simple fringe projection profilometry (FPP) system for the analysis of semi-fossilized samples.

## 2. Description of some phase shifting algorithms (PSA)

Many authors have developed a variety of methods for PSA, such as averaging with existing algorithms and solving for the roots of a characteristic polynomial [7], data-sampling windows, Fourier analysis, least-squares, etc. The use of the TPS is based on the addition of a careful phase change that is projected to the target surface. It also involves analyzing data from each pixel independently of all other pixels in the frame. This technique is based on the work of Carré as shown in Eq. (1), considering that the problem of phase-shifter miscalibration is dealing with treating the phase shift as one more unknown variable.

## 3. Designing an eight-frame phase shifting algorithm

In phase shifting interferometry the ideal intensity $I\left(x,y,t\right)$ for k = 1, 2, 3… *M* of each interferogram recorded by a CCD detector can be expressed as [6–8],

*x*and

*y*denote the pixel position; $a\left(x,y\right)$is the background illumination; $b\left(x,y\right)$ is the contrast of interference fringes (amplitude), and $\phi \left(x,y\right)$ represents the phase. Meanwhile, the temporal carrier ${\omega}_{0}$ is a linear phase shift among this set of interferograms and it is determined by how fast the phase reference wave is changing.

By modeling the fringes with the equation shown in [5], the PSP problem is usually reduced to four steps:

- I.
*M*images are captured with several phase shifting among them. - II. To choose or design a specific M-frames phase shifting algorithm (PSA) to process the set of M images to obtain the wrapped phase.
- III. An unwrapping algorithm to recover the desired phase is designed.
- IV. A texture is applied to the obtained phase to exhibit the desired target.

It should be highlighted that an eight-frame algorithm corresponds to a filter with seven independent parameters to recover the desired phase and to compensate for some errors. Then, at least two of those parameters are necessary to eliminate the D.C. component and the fundamental frequency. Two additional conditions are used to compensate the linear bias variation and the linear phase shift detuning error [16]. The remaining three conditions are used to compensate other errors generated by other effects, like a non-linearity response and to obtain a better SNR. After testing well-known filters of three, four, five, six, seven and more steps, we observed that using seven or more frames the resolution obtained was in the parameters required for this type of targets. The response did not improve substantially when using nine or more frames. Yet, the decisive point was that when analyzing the histogram of the estimated frequency of tuning this was broadband and tuned on $\pi /4$ where the filter of Hibino had the best response closer to the one expected. However, this filter was not designed for the conditions of our experiment [17], thus for that reason we generated a robust filter range, designed specifically for the needs of our experiment: a broadband filter, tuned on $\pi /4$, insensible to bias variations and detuning errors, specially linear and a good SNR.

The estimated phase of any quadrature filter with an M^{th} order is given by [6,7,16,18],

*N*and

*D*are the desired numerator and denominator row vectors. In a previous work [1,16,18], it was proved that the Fourier transform $H\left(\omega \right)$ of this filter is,

*M-1*frequencies that are the necessary conditions to be a specific filter [16,18]. Therefore, from Eq. (6), the general case for

*M*= 8, the corresponding eight-frame algorithm is,

*N/D*with symmetric coefficients is obtained from the expression [16,18]

*π/4,*it is well known that a filter that eliminates harmonics corresponds to the cut-off frequencies ${\alpha}_{k}$ =

*0, π/4, π/2, 3π/4, π, 5π/4*and

*3π/2*[1,11,16,18]. That is, to be a quadrature filter, the two necessary cut off frequencies are zero and the phase step π/4 [1]. However, in a least-squares procedure the harmonics are also eliminated. Therefore, from Eq. (10) and after solving some algebra operations the obtained eight-frame algorithm is,

*M = 8*, the result equivalent to Eq. (10) becomes,

*α*as the ideal phase step to tune the filter, then a new (

*n + m*-1) frame filter is obtained from two individual filters as shown in Eq. (7). As mentioned before, the design of a tunable filter allows this case to be extended further to an eight-frame filter, which allows the selection of the data to be removed. Furthermore, we select polynomial roots implying that the filter must suppress frequencies in Fig. 1 . To assure that the filter eliminates harmonics, undesirable frequencies and the systematic errors involved, we first propose a filter to deal with harmonics, mainly as in Eq. (13). Also, we define a filter that is able to handle an optimal SNR and linear detuning errors as

*0, π/4, π/4, π/2, 3π/4, 3π/4, π*(from now known as mainly detuning error filter or MDE filter). Applying the same method to design a filter that compensates mainly bias errors (MBE filter) associated to the system, we calculate a filter considering cut off frequencies in

*0, 0, π/4, π/2, 3π/4, π, π,*and finally a filter centered in

*0, 0, π/4, π/4, π/2, 3π/4, π*to compensate detuning and bias errors (DBE filter). Graphic representation of these filters with their cut off frequencies is shown in Fig. 1 according to [6].

Then, applying the above procedure for the rest of the filters shown in Fig. 1, from Eq. (10) the filter shown in Fig. 1(b) is,

Notice that the obtained filters above satisfied a wide range of requirements for any experiment. That is, considering that the ideal quadrature filter is a step function, the ratio of the area under the curve between the right and the left side of the graphic, starting from zero, gives us an idea of how the filter works, and why the best approximation for this application is depicted as the green line corresponding to a PSA compensating detuning or MDE. Additionally, we can observe the response for the phase shift detuning, bias, and harmonics (Fig. 3 )

This graph verifies the design of the filters MBE and DBE, having a double cutoff frequency at zero and therefore it is sufficient to compensate the bias error. In the same way, the green trace (or MDE filter) in the tuning frequency $\pi /4$ has a similar response and this makes it insensitive to linear detuning error. Furthermore, linearity is given by the symmetry of the function in $\pi /2$ and the filter with the lowest area under the curve has better SNR. Since each filter is defined to counteract a specific error, we cannot choose one of the algorithms as the best for a given application. Indeed, all proposed filters are good and efficient to compensate for the error by which they were designed. The phase shift detuning error is depicted in Fig. 4 .

As expected, when the difference between the ideal and the observed carrier is zero, the phase error is zero. Having a carrier variation between −0.2 and 1.6 rad, the estimated phase error is less than 0.03 rad. This implies that the filters in Figs. 1(b)–1(d) are broadband filters, where the best performance is from the filters MDE (Fig. 1(b)) and BDE (Fig. 1(d)).

## 4. Evaluation of the algorithms

#### 4.1 Simulation

In order to evaluate the robustness of the proposed filters, we use a simulation of a profile to quantify the error of those filters and to compare such results with a single four-frame filter and a Hibino filter (such comparisons were made among Surrel [11], Hibino [13] and de Groot [19] filters; however, the Hibino filter shows a better curve fitting with the parameters included in these simulations). First, we simulate a profile by using the Eq. (5). Random noise is then added as detuning error in the form of $\omega =\omega +\Delta \omega $and bias error as $b=b+\Delta b$ to the fringe pattern. Increasing the detuning error until a 10% and the bias error until a 20%, we started to see the different trends of each of the filters as shown in Fig. 5 .

The best filter observed, the MDE filter, has an error average of 0.2% and an average of deviation of errors of 0.4%. Table 1 shows the average of errors and deviations of this simulation.

From Table 1 we can notice that the function that has a lower error rate is, as expected, the least-squares filter, this is because the function is designed to precisely obtain the least square error of the simulated points. Nonetheless, the MDE filter has better performance with less error because despite it has a slightly higher average error than the least-squares filter, in the deviation of data we can see that it is much higher for the latter. In the same manner, now we simulate a more complex function, as the well-known peak function as shown in Fig. 6 . In this case we use a 15% of bias and 10% of detuning randomly added to the original plot.

Then, to obtain the standard deviation and average error indicators, first we subtract the graphic obtained with random errors from the ideal graphic without errors, followed by an average of all the points to finally obtain the average error of that filter in the simulation. For deviation, initially we have to get the standard deviation of all the points to then average them. Table 2 shows the results obtained. The closer is the correlation coefficient to one, the more it resembles to the original function.

From Table 2, it is observed that the lesser average error corresponds to the MDE filter; however, a four-step filter was better than MBE filter due to the fact that while the more steps are acquired, the error introduced in the algorithm for obtaining the wrapped phase also increases. From the simulations we can conclude that for cases where detuning errors or bias errors are present, proposed tuned filters on $\pi /4$ give us a very good method to retrieve accurate information of an object despite the experimental random noises. Finally, from the Pearson’s coefficient we obtain the correlation between the original sample and the estimated phase and as we expect the filter MDE was the best choice again.

#### 4.2 Experiment

A series of eight frame experimental patterns were obtained from a fringe projection profilometry system [20,21] which were used to test several eight-frames filters and the algorithms implemented (the system is not discussed in this paper). To evaluate the algorithms in the experiment and identify the most suitable, we used a semi-fossilized sample of a rodent mandible with a length of 20 mm, provided by the Instituto Nacional de Antropología e Historia (INAH). Semi-fossilized materials were selected because of their intrinsic challenges [5], such as the need of a non-contact experiment due to its fragility and risk of modern DNA contamination associated with sample handling, its extremely small size, and the need to preserve the ancient DNA intact. A set of eight images of this sample were digitalized and selected considering the level of quality required (Fig. 7 ).

The intensity patterns were recorded by using a CCD camera with a resolution of 1280x1024 pixels in grayscale and a National Instruments data acquisition card. The measurement process was controlled by NI LabVIEW, and the system was calibrated in gamma and sinusoidal projection through this software. The images were processed with MATLAB software applying several algorithms; however, systematic errors and miscalibrations were still present. Since each displacement in phase has a partial error shift when the relative phase increment occurs, it is necessary to estimate the phase shifts error.

## 5. Results and discussion

By using Carré’s phase estimator, the error in the phase shift among interferograms was revealed. Therefore, the need to design a set of eight-frame algorithms for phase shifting profilometry (PSP) immune to systematic errors, such as higher harmonics in the intensity signal, detuning, nonlinear response of photodetectors, random intensity noise and vibrations [17] in these measures becomes essential as shown in Fig. 8 .

For the development of the best mathematical method for processing the phase shift, several comparisons have been made considering the immunity to errors and miscalibrations of each of the algorithms mentioned above (Fig. 9 ). Accordingly, the width and shape of the histogram represent the variation of phase shifts.

For these interferograms, it is evident that Carré’s approach is susceptible to higher-order harmonics in the signal [9]. Filters reported previously in [12] have a good fitting to phase miscalibrations, but in the case of turbulences this method is not recommendable. Results from well-known algorithms, as proposed by Schmit (Fig. 9(a)) and Carré (Fig. 9(b)), show a high sensitivity to systematic errors. Hibino and Surrel (Figs. 9(c)–9(d) respectively) methods are more sophisticated in correcting errors, but there are some inconsistencies in the form of the wrapped phase. Least-square filter in Fig. 10(a) is enough for most applications, because many do not require a maximum level of immunity to detuning and miscalibrations of the phase shifts. Methods proposed have the robustness to manage interferograms with physical limitations and the intrinsically experimental errors present in this kind of measures. The first one is the well-known least-square eight-frame filter, the second corresponds with an algorithm specialized in compensating the phase shift linear detuning error, while the third algorithm is designed to compensate the linear bias variation error and the fourth is an algorithm that simultaneously compensates both kind of errors. However, a successful result depends on two main facts, the algorithm and the kind of errors and noise present in the set of frames measured.

In addition, the evaluation of intensities can be made by plotting a linear slice across interferograms and analyzing the profile of the phase as in Fig. 11
. SNR evaluation of our proposed methods is significantly lower. The MDE algorithm shown in Fig. 10(b) with cut off frequencies in *0, π/4, π/4, π/2, 3π/4, 3π/4, π* has the better immunity to experimental errors in this FPP arrangement and these samples in particular.

These results are consistent with the unwrapped heights map, and it gives the right step for the analyses. A textured map of unwrapped phase is shown in Fig. 12 . When the algorithm is limited in suppressing errors, we generally need more expensive equipment for the experiment to compensate them. On the other hand, having a good algorithm and considering the particular experimental conditions, more economic equipment will be enough to cover the needs of precision and accuracy of the experiment.

Most eight step algorithms are centered at π/2 [22]; however, as the error distribution is greater for each step, a filter centered in a further step (π/4 in this case) and immune to detuning (given the dispersion of the phase steps in Fig. 8) is needed to compensate experimental errors. In our experiment, objects under study were small (~20mm), and the technique proved to be a simple, cheap, and flexible method to be applied in areas that use phase shifting interferometry, but also other fields like paleontology, anthropology, and the sort. For instance, studies that need high resolution images for a posteriori scientific analyses, as well as morphological studies that involve ancient DNA preservation and reconstruction of semi-fossilized samples at different scales, are some of the applications and advantages of the described technique.

## 6. Conclusions

We have obtained four novel filters corresponding with eight-frame algorithms to solve the problem of the reconstruction of three-dimensional surfaces in semi-fossilized materials. They are able to estimate the phase shifts by minimizing the high-order correlations between the reconstructed phases and the reconstructed amplitudes of the reference wave, which are introduced by phase shift errors in PSP. The proposed algorithms are not reported previously, they were tested experimentally besides simulations. Too they were compared with traditional eight step PSP methods. In particular, the MDE filter reduces significantly the influence of the detuning error and harmonics. This algorithm has several advantages over similar ones. First, it allows us to analyze materials with several optical properties involved in miscalibrations and phase errors. Second, it is suitable for studies that involve FPP systems. The effectiveness and improvement of the proposed algorithms and procedures have been supported and confirmed by our simulations and experimental results.

## Acknowledgments

This work was supported by the Mexican National Science and Technology Council (Consejo Nacional de Ciencia y Tecnologia, CONACYT) through scholarship grants 175434 and 207796, project #180449, and CONACYT-MAE Program 146523. The authors thank the Instituto Nacional de Antropología e Historia (INAH) for the loan of the semi-fossil sample (C.A.401-36/0360). We also acknowledge the help granted by MDD J. J. Lozano in the revision of this paper.

## References and links

**1. **J. F. Mosiño, D. M. Doblado, and D. M. Hernández, “A method to design tunable quadrature filters in phase shifting interferometry,” Opt. Express **17**(18), 15772–15777 (2009). [CrossRef] [PubMed]

**2. **K. Larkin, “A self-calibrating phase-shifting algorithm based on the natural demodulation of two-dimensional fringe patterns,” Opt. Express **9**(5), 236–253 (2001). [CrossRef] [PubMed]

**3. **P. D. Ruiz, J. M. Huntley, and G. H. Kaufmann, “Adaptive phase-shifting algorithm for temporal phase evaluation,” J. Opt. Soc. Am. A **20**(2), 325–332 (2003). [CrossRef] [PubMed]

**4. **H. Katterwe, “Modern approaches for the examination of toolmarks and other surface marks,” Forensic Sci. Rev. **8**, 45–72 (1996).

**5. **T. M. Kaiser and H. Katterwe, “The application of 3D-microprofilometry as a tool in the surface diagnosis of fossil and sub-fossil vertebrate hard tissue. An example from the pliocene upper laetolil beds, Tanzania,” Int. J. Osteoarchaeol. **11**(5), 350–356 (2001). [CrossRef]

**6. **Y. Surrel, “Design of algorithms for phase measurements by the use of phase stepping,” Appl. Opt. **35**(1), 51–60 (1996). [CrossRef] [PubMed]

**7. **D. W. Phillion, “General methods for generating phase-shifting interferometry algorithms,” Appl. Opt. **36**(31), 8098–8115 (1997). [CrossRef] [PubMed]

**8. **J. Schmit and K. Creath, “Extended averaging technique for derivation of error-compensating algorithms in phase-shifting interferometry,” Appl. Opt. **34**(19), 3610–3619 (1995). [CrossRef] [PubMed]

**9. **K. Creath, “Temporal phase measurement methods,” in *Interferogram Analysis*, D. W. Robinson and G. T. Reid, eds. (Institute of Physics, 1993).

**10. **J. M. Huntley, “Automated Analysis of Speckle Interferograms,” in *Digital Speckle Pattern Interferometry and Related Techniques*, P. K. Rastogi, ed. (Wiley, 2001).

**11. **Y. Surrel, “Phase stepping: a new self-calibrating algorithm,” Appl. Opt. **32**(19), 3598–3600 (1993). [CrossRef] [PubMed]

**12. **K. Creath and J. Schmit, “N-point spatial phase measurement techniques for nondestructive testing,” Opt. Lasers Eng. **24**(5-6), 365–379 (1996). [CrossRef]

**13. **K. Hibino, B. F. Oreb, D. I. Farrant, and K. G. Larkin, “Phase-shifting algorithms for nonlinear and spatially nonuniform phase shifts,” J. Opt. Soc. Am. A **14**(4), 918–930 (1997). [CrossRef]

**14. **C. Rathjen, “Statistical properties of phase-shift algorithms,” J. Opt. Soc. Am. A **12**(9), 1997–2008 (1995). [CrossRef]

**15. **C. S. Guo, L. Zhang, H. T. Wang, J. Liao, and Y. Y. Zhu, “Phase-shifting error and its elimination in phase-shifting digital holography,” Opt. Lett. **27**(19), 1687–1689 (2002). [CrossRef] [PubMed]

**16. **J. F. Mosiño, J. C. Gutiérrez-García, T. A. Gutiérrez-García, and J. M. Macías-Preza, “Two-frame algorithm to design quadrature filters in phase shifting interferometry,” Opt. Express **18**(24), 24405–24411 (2010). [CrossRef] [PubMed]

**17. **P. D. Ruiz, J. M. Huntley, and G. H. Kaufmann, “Adaptive phase-shifting algorithm for temporal phase evaluation,” J. Opt. Soc. Am. A **20**(2), 325–332 (2003). [CrossRef] [PubMed]

**18. **J. F. Mosiño, J. C. Gutiérrez-García, T. A. Gutiérrez-García, F. Castillo, M. A. García-González, and V. A. Gutiérrez-García, “Algorithm for phase extraction from a set of interferograms with arbitrary phase shifts,” Opt. Express **19**(6), 4908–4923 (2011). [CrossRef] [PubMed]

**19. **P. Groot, “Derivation of algorithms for phase-shifting interferometry using the concept of a data-sampling window,” Appl. Opt. **34**(22), 4723–4730 (1995). [CrossRef] [PubMed]

**20. **S. S. Gorthi and P. Rastogi, “Fringe projection techniques: whither we are?” Opt. Lasers Eng. **48**(2), 133–140 (2010). [CrossRef]

**21. **S. Ma, C. Quan, R. Zhu, and C. J. Tay, “Investigation of phase error correction for digital sinusoidal phase-shifting fringe projection profilometry,” Opt. Lasers Eng. **50**(8), 1107–1118 (2012). [CrossRef]

**22. **J. A. N. Buytaert and J. J. J. Dirckx, “Study of the performance of 84 phase-shifting algorithms for interferometry,” J. Opt. **40**(3), 114–131 (2011). [CrossRef]