1. Introduction

Interferometry has become a powerful tool in metrology owing to the sensitivity of optical phase to a range of relevant measurands [1, 2]. In single frequency interferometry the fractional fringe value can be calculated by either phase stepping or Fourier transform techniques [3], but the intrinsic periodic nature of the fringe function means that fringe order information is lost. This fundamental limitation must be resolved in order to obtain absolute phase measurement. Multi-frequency interferometry has provided a solution for fringe order identification with classical analysis performed using the method of excess fractions [4, 5, 6, 7]. In this technique the fractional fringe values measured at each frequency are used to form a set of simultaneous equations that may be solved for the integral fringe order. The measurements are absolute within an unambiguous measurement range (UMR) that is maximized by ensuring that none of the measurement wavelengths share common factors, i.e. they are pair-wise relatively prime. This gives a single coincidence of wavelength at all measurement frequencies over the UMR [7].

Alternative approaches for absolute phase measurement have been proposed, based on multi-wavelength interferometry with beat frequency processing to generate a synthetic wavelength from a coincidence of phase [8, 9, 10]. It was found that the UMR of these techniques was limited by phase measurement noise. Techniques to further extend the UMR beyond the beat wavelength have also been reported and are relevant when the optical measurement wavelengths are well separated [10]. Recently, we introduced a robust optimization process for frequency selection in multi-wavelength interferometry where a geometric series of synthetic wavelengths is defined to maximize the overall process reliability [11]. This imposes a different frequency selection criterion compared to the method of excess fractions and hence both approaches will find particular applications depending on specific experimental constraints and the available optical sources.

The processing time to find the absolute phase using the method of excess fractions is linearly dependent on the UMR and becomes computationally prohibitive for data from full-field interferometers. Important applications include transient events and real time control of dynamic forming [12]. Processing times may run to several days for typical data sets [7], hence faster processing algorithms are desired.

A solution to this computational problem has been proposed based on the Chinese Remainder Theorem (CRT) from number theory [13, 14, 15, 16] where wavelength selection is based on pair-wise relatively prime wavelengths. Unlike the excess fractions method a direct calculation of the measurand is produced using modulo arithmetic. To apply the CRT method, wrapped phase values are scaled to integers and multiplied by a large integer wavelength product, hence the approach is very susceptible to phase measurement noise [14, 15]. Agurok [14] demonstrated that considering a reduced UMR could accommodate small errors in the wrapped phase measurements. Alternatively, Takeda et al. [15] proposed a method where phase offsets were used to move data away from integer boundaries and it was recognized that fringe order errors form a structured distribution. Spatial comparisons between neighboring pixels were checked against this distribution and if a match was found, the values were corrected. This error correction procedure assumed relatively clean data, i.e., where the integer wrapped phase errors have a maximum magnitude of 1. Zhong et al. also reported a modified form of this approach based on an irrational ratio of wavelengths [16]. In all cases data was presented for a two-wavelength system over a limited number, <10, of fringe orders.

In this paper we describe a new modified form of the CRT algorithm to accommodate phase noise in the input data whilst maintaining a low computational burden. The algorithm is applicable to a generic r-wavelength interferometer and is based on a computationally efficient search over an (r-1) dimensional matrix. The performance of the modified CRT algorithm has been compared to an excess fractions solver by evaluating the reliability of the measurements obtained and the process time. Experimental results are presented from a 3-wavelength full-field shape measurement system where the data has been processed using the modified CRT algorithm.

2. Theory

Consider an interferometer using r measurement wavelengths given by λ ₁, λ ₂,…,λ_r . At each wavelength a wrapped phase measurement is obtained, φ_i , i=1,…,r. To apply the CRT, the wavelengths must be expressed as integers by applying an appropriate scaling factor. Similarly, the wrapped phase values must be cast as integers by multiplying by λ_i /2π and applying the same scaling factor. The CRT states that given the system of congruences [13, 16]:

there is a unique solution for x, the absolute distance required, modulo Λ=λ ₁ λ ₂…λ_r provided that the greatest common divisor of all pairs of wavelengths is 1, i.e. gcd(λ_i , λ_j )=1∀i≠j. The integer x is a simultaneous solution of the r congruence’s given by:

where Λ _i =Λ/λ_i and Λ′ _i is the inverse of Λ _i modulo λ_i . The Λ′ _i terms are calculated once for any set of measurement wavelengths using the extended Euclidean algorithm [13]. Therefore, the measurand is found by direct calculation and the processing time is independent of range.

The sensitivity of the CRT algorithm to phase measurement noise in the φ_i values can be illustrated by means of an example. For a 3-wavelength system with measurements at 473, 532, 633 nm and considering 1 nm as the integer unit, the Λ _i Λ′ _i terms are -23909676, -53594211 and 7750388 respectively. Hence a 1 nm error in φ_i , corresponding to a phase resolution of <1/500^th of a fringe, causes a very large change in the calculated value of x.

2.1 Modified CRT Algorithm

In a practical interferometer the measured phase $\widehat{\phi }$_i, contains noise, Δφ_i , such that $\widehat{\phi }$_i=φ_i +Δφ_i . This leads to an uncertainty in the measurand, Δx, given by:

where the noise present in each wrapped phase measurement is expressed as an integer. A matrix of values of Δx can be generated around the measured value corresponding to $\widehat{\phi }$_i. The noise distribution at each measurement wavelength will have an upper and lower bound; the actual form of the distribution is immaterial to the operation of the algorithm. By defining n as the maximum absolute value of the upper and lower bounds for each wavelength we can define an interval with (2n+1) elements containing these bounds. The r-dimensional matrix of Δx values is then setup such that Δφ_i =-n,-n+1,…,0,…,n-1, n along each axis, with the total number of elements in the matrix given by (2n+1) ^r . By adding the Δx values to the initial calculation for x corresponding to $\widehat{\phi }$_i, a matrix of possible solutions for the measurand is determined. Each element in the matrix is then potentially the correct solution for the value of x that would be obtained in the noise free case. By searching the matrix and using a sufficiently restricted UMR it is possible to find the correct solution. The minimum separation between pairs of values within the matrix then defines the reduced UMR of the sensor.

There is a compromise between the probability of finding the correct solution and the associated reduction in UMR. As the search space increases the range of values taken by the sum of integer multiples of Λ _i Λ′ _i also increases, giving a greater possibility that the elements of the matrix contain similar values. Hence the UMR is reduced. However, as the size of the search space is increased there is a greater probability that the correct solution is contained within the search space, assuming that the phase noise distribution is non-uniform. If there is a wide choice of measurement wavelengths available then careful selection allows the reduced UMR to be optimized.

In implementing the modified CRT algorithm the emphasis has been placed on calculating the correct fringe order, i.e., NINT(x/λ_i ), where NINT denotes taking the nearest integer. It has been observed when examining an r-dimensional matrix, corresponding to r measurement wavelengths that all possible fringe order solutions exist in a matrix with (r-1) dimensions. This occurs because the sum ${\displaystyle \sum _{i=1}^r}$Λ _i Λ′ _i =±1(modΛ). Therefore the computational efficiency of the process can be improved by searching the reduced (r-1) space. To ensure that the same set of solutions is obtained over (r-1) dimensions, the span of the search must be increased giving a total search space of order (4n+1)^r-1. Overall there is a significant computational benefit for low numbers of measurement wavelengths, this gain increases as n increases. For example, in a 3-wavelength system, the search space decreases by a factor of 1.5 for n=2, and by a factor of 4 for n=7. A summary of the algorithm is given in appendix A.

3. Simulations

The performance of the modified CRT algorithm has been evaluated using a computational model with 3 measurement wavelengths at 517, 558 and 641 nm. For each wavelength, the model generates a linear object phase profile containing 10000 data samples from which wrapped phase distributions are calculated. A Gaussian noise model is appropriate for most interferometers where there is sufficient illumination intensity [17] and furthermore defines the probability that the correct solution will be found within a given search matrix size. In this model, noise is added to the phase distributions with zero mean and a standard deviation of either 2π/200 or 2π/500 radians. The resultant phase distributions are then wrapped into the interval -π to +π to represent the measurements obtained from a real sensor. The wrapped phase values are scaled to integers where 1 unit represents 1 nm prior to calculating the absolute range from equation 2. The span of the search space was selected to encompass at least ±2 standard deviations of the added phase noise in order to ensure that the search space contains the correct solution at >95% of the data samples. The results of the simulations are summarized in Figs. 1 and 2 for 2π/200 and 2π/500 phase noise respectively. Each figure contains data from 3 algorithms: a conventional excess fractions solver, an implementation of the standard CRT algorithm and the modified CRT algorithm. For each algorithm two data sets are plotted; the time taken (running on an 800MHz processor in Matlab) and the percentage of data samples giving the correct fringe order, a measure of the reliability of the algorithm. The dependence of these parameters on the depth range of the object profile is shown in Figs. 1 and 2. It should be noted that the theoretical UMR for these wavelengths is approximately 4 times larger than the maximum object range examined here.

Figures 1 and 2 show that the time taken to calculate the solution for the excess fractions solver has a linear dependence with range. This was expected as the object range simulated provides a priori knowledge of the depth range to search over for integer fringe order solutions at each measurement wavelength. The small deviations from the linear trend over small object ranges are due to increased processing overheads in loop implementation within the code. For excess fractions the results show a decrease in the percentage of data samples giving the correct fringe order as the range is increased. A reliability of >90% was achieved for a range of up to 250 microns with a phase noise of 2π/200 radians and for a range of up to 2000 microns for a phase noise of 2π/500 radians. The conventional CRT algorithm in both cases shows an unacceptable reliability of a few percent. As expected the time taken by the CRT algorithm is constant with respect to range.

For a phase noise of 2π/500 radians the modified CRT algorithm shows a processing time advantage compared to the excess fractions solver for ranges greater than 20 microns. The modified CRT is >420 times faster than the excess fractions solver at a range of 2000 microns where both algorithms exhibit similar reliabilities in excess of 98%. For an object range >2000 microns both methods give reliabilities of <90% that is unlikely to be acceptable in practical applications. For a phase noise of 2π/200 radians the processing time advantages are reduced with a time benefit of a factor of 10 at a range of 250 microns. Above this range the process reliabilities for both algorithms fall to unacceptable levels.

 

Fig. 1. Comparison between Excess Fractions, CRT and Modified CRT algorithms with phase measurement noise of 2π/200 radians.

Download Full Size | PDF

 

Fig. 2. Comparison between Excess Fractions, CRT and Modified CRT algorithms with phase measurement noise of 2π/500 radians.

Download Full Size | PDF

4. Experimental results

The modified CRT algorithm has been validated using experimental data from a practical full-field fringe projection system, shown schematically in Fig. 3. A commercial data projector was utilized to generate the fringe patterns. It was found that a non-linear gray scale response existed between the intensities transmitted by the data projector and those received on a 12-bit digital CCD camera. A cubic polynomial function was applied to the intensities transmitted by the projector to compensate for this non-linearity. The measurement phase noise was evaluated using a flat test object and found to follow Gaussian statistics with a standard deviation on of ~2π/250 radians.

To apply the modified CRT algorithm using the full-field data projector, the measurement wavelengths were synthesized to the wavelengths in pixels of the projected fringes. For this evaluation, wavelengths of 14.4, 19.1 and 23.3 were selected so that when multiplied by 10 they are pair-wise relatively prime. The search space along each axis was set such that approximately ±3 standard deviations of the measured phase noise would be accommodated, i.e. a reliability of >99.73% will be obtained through the search process in the modified CRT algorithm. Initial measurements were made on a flat plate and the results validated against the expected fringe order distribution at the shortest wavelength. A reliability of 100% was obtained over the 71 fringe orders in the image.

 

Fig. 3. Experimental arrangement.

Download Full Size | PDF

Measurements from a typical test object, in this case a rapid prototype model of an engine port and combustion chamber, is presented. The wrapped phase distribution, at the shortest wavelength, and the corresponding unwrapped distribution are shown in Fig. 4(a) and (b) respectively. A low modulation check has been made in each of the wrapped phase maps to identify pixels where there is insufficient light scattered by the object. These points are then excluded from analysis by the modified CRT algorithm. As expected, the absolute phase has been determined at each pixel independently.

 

Fig. 4. (a) Wrapped phase map from an engine port model.

Download Full Size | PDF

 

Fig. 4. (b) Unwrapped phase map from an engine port model.

Download Full Size | PDF

5. Summary

In this paper we have introduced a modified form of the Chinese remainder theorem algorithm for processing multi-frequency interferometric data that is both robust and time efficient. The algorithm is applicable to a generic interferometer using r measurement wavelengths. A search scheme to compensate for phase measurement noise has been implemented where the search space has been reduced to (r-1) dimensions. The computational benefit of the process increases as the search space is increased. It has been demonstrated that the modified CRT algorithm applied to a 3 wavelength system reduces the processing time by a factor of >420 compared to the method of excess fractions (Fig. 2) whilst achieving similar process reliability. Therefore, in the future this technique offers the potential to make processing times practicable for full-field shape measurement from high-speed transient events or in real time process control applications. Furthermore, for a phase noise of 2π/500 radians the simulations with both the modified CRT and excess fractions solvers show that only ~1% of the theoretical UMR can be measured with >98% probability of obtaining the correct measurements.