Application of matched digital filters to noisy fringe-patterns from complex wavefronts

Application of matched digital filters to noisy fringe-patterns from complex wavefronts J.Caum, J.Arasa, S.Royo and M.Ares Centre for Sensor, Instrumentation and Systems Development, Technical University of Catalunya (UPC-CD6) Rambla Sant Nebridi 10 E08222 Terrassa Spain ABSTRACT A new technique for processing of noisy fringe patterns with complex fringe shapes is presented . The technique is based in exploiting the capabilities of digital filters, treating the fringe-pattern as a signal to be processed. The main parameters required for designing the filter, and the application of the design methodology to a particular experimental fringe pattern are described in detail. To conclude, the designed filter is applied to experimental noisy fringe-patterns both from complex and simple wavefronts, which are processed without a single change in the parameters of the algorithms. The procedure is completel y general, with a very small computational cost, and may b e extended to any nois y fringe-pattern to be processed given the fringes are still visible. Keywords: Fringe processing, metrology, optical testing, discrete digital filters INTRODUCTION Many optical shop techniques rely on the particular fringe processing technique applied, so fringe processing has been a widely explored topic where much research has been carried out. Fringe patterns have become a valuable source of information on different techniques, from interferometry to holography, yielding its processing a key part of the complete metrolog y strateg y o f the tool designed 1 . In particular, the location of the centre of the fringe, and the possible change in phase introduced in the fringe pattern by the image processing algorithm, among other factors, may noticeably disturb the final result of the measurement. In one of the usual approaches, when only the position of the maxima and minima of the fringes are of interest, they are usually extracted through classical morphological procedures to yield a final skeletonized pattern where a single line sums up the information contained in the whole fringe. Complete phase information on the signal may then be extracted from the analysed fringe-pattern. The classical approach, described elsewhere 2 , involves a binarization procedure followed by a skeletonization procedure based in eroding the binarized fringe pattern; Fourier-transform based methods have also be en proposed. H o w e ver, these methods have three main drawbacks: they are quite intensive computationally, espe cially for the eroding step in complex wav efront shapes, as they are morphological procedures applied to the whole fringe pattern; they are dependent on a bundle of parameters which change from fringe-pattern to fringe-pattern, as for instance the binarization level, or the parameters determining the skeletonization strategy; but, over all of them, they rely on a binarization across the whole field, which is a process highly dependent on intensity changes of a given fringe or on the image intensity bias across the field of view. This often leads to fringe discontinuities which cause errors in the interpretation of the results. Although some alternative approaches for complex fringe shapes have been proposed 3 , noisy and complex-shaped fringe-patterns still defy automatic image processing, as demonstrates the important amount of literature already devoted to the topic 4 . The method we present overcomes all these drawbacks, by applying a digital filter matched to the spatial frequency of the fringe-pattern to each of the lines in the image. B y u sing a loworder IIR-ARMA filter a very g o od attenuation with negligible phase change is introduced, allowing a same filtering strategy to be applied to simple and complex fringe patterns without any parameter adjustment but the cut-off frequenc y desired f or the filter. With the remaining data, both the detection of the maxima of the fringes or the phase measurement may be obtained from the image. The designed matched filter is applied to experimentally measured wavefronts using a digital Ronchi deflectometry s etup 5 and applied Please verify that (1) all pages are present, (2) all figures are acceptable, (3) all fonts and special characters are correct, and (4) all text and figures fit within the margin lines shown on this review document. Return to your MySPIE ToDo list and approve or disapprove this submission. 5856-83 V. 1 (p.1 of 11) / Color: No / Format: A4 / Date: 2005-03-29 13:19:04 SPIE USE: ____ DB Check, ____ Prod Check, Notes: to the positioni ng of the maxim a of the fri nges, without loss of gener al ity regard i ng app lica tio n to any other experi me nt al technique or to phase extraction. With this purpose, the paper describes the basic c oncepts and nomenclature regarding digital filter design along Section 2, as a short introduction to the nomenclature of digital signal processing may be of much help to non-electronic researchers developing optical metrology applications. The complete design procedure of a given matched filter is described along Section 3. Finally, the application of the filter to both simple (spherical) and complex (progressive power) wavefronts is presented in section 4, showing the validity of the approach and its general application both to simple and complex fringe patterns without parameter adjustments. Section 5 briefly describes the conclussions drawn out. 2. INTRODUCTION TO DIGITAL FILTER DESIGN Digital filter design is usually de v oted to the proces sing of time-dependent periodic signals. In our case, a fringe-pattern will become a periodic signal with an spatial frequency determined by the size of the CCD array and the number of fringes present in the image. Our purpose is to design a digital filter in the spatial domain from the well-stablished set of procedures involved in digital filter design. This pushes to a short introduction on the theory of digital filter design, where only the more relevan t aspects from the poin t of view of fri nge-pa tte r n p rocessing will be presented, as fa r as the re is a wide variety of high-quality literature on the topic 5 . The main goal of our approach is to reduce the high (spatial) frequency components of a signal. The signal to be processed is each of the individual pixel rows in the image, which will be processed in sequence until the complete fringe pattern has been analys ed. Obtaining a smooth signal with low-frequenc y information obvi ousl y pushes to a lo wpass digital filter, where the attenuation of the signal under the cut-of f frequenc y Ω p is small (typicall y under 3dB), and the attenuation after the minimum attenuated frequenc y Ω a is very important. The slope in the transition region and the maximum attenuation attained depend on the value of the order of the filter N . A typical representation of the tranfer function H(j Ω ) is presented in fig.1. Fig.1: Typical transfer function for a general low-pass filter: Ω p is the cut-off frequency, Ω a is the minimum attenuated frequency, δ p is the amplitude of the modulus of the transfer function in the passing frequencies and δ a is the s ame va lue in th e at ten ua ted frequencies Please verify that (1) all pages are present, (2) all figures are acceptable, (3) all fonts and special characters are correct, and (4) all text and figures fit within the margin lines shown on this review document. Return to your MySPIE ToDo list and approve or disapprove this submission. 5856-83 V. 1 (p.2 of 11) / Color: No / Format: A4 / Date: 2005-03-29 13:19:04 SPIE USE: ____ DB Check, ____ Prod Check, Notes: This type of filters may be implemented through different approaches, being filters from the Butterworth, Chebyshev and the elliptical families the most usual. In our case, we have selected the Butterworth approximation for filter implementation, which has a transfer function H(j Ω ) with a modulus usually expressed as () N p c jj jH 2 2 /1 1 )( ΩΩ+ =Ω (1) where j stands for the imaginary constant, Ω stands for frequency and Ω p is the abovementioned cut-off frequency of the filter. The Butterworth approach has been selected because it does not present undulations in the modulus of the transfer function both in the passing and in the attenuated frequency ranges (this is to say, it is monotonic over both frequency ranges). For fringe pattern analysis, unit gain in the passing range of frequencies is convenient in order to keep the contrast of the signal constant. Another main issue which pushed us to select this family of filters is the change it introduces in the phase of the filtered signal, which is always smooth and continuous. An additional aspect to be considered is the order of the filter N, as the higher the order of the filter, the steeper the phase change occurring around the cut-of f frequenc y, but also the higher the computational co st of the filtering. A loworder filter with the capability o f properly filtering the signal is the obviously the optimum solution. One of the main issues to be considered is the change in phase introduced in the signal, as it appears to be a limiting factor which, after the application of the filter, will displace the positions of the maxima of the fringe in a given di rection . T his disp lacemen t is rel ated i n dig it al filte r desig n to t he g roup del ay τ ( Ω ) introduced by the filter, defined as Ω Ω −=Ω d jHd c ))((arg )( τ (2) Up to the moment, the continuous frequency approach has been asumed. However, our signal will be spatially discrete, so we will leave the continuous domain to enter the discrete frequency domain. Obviously, our next step is to define which type of transfer function is more adequate for the fringe-processing problem. A simple classification divides digital filters into two families with radically different transfer functions: a) Finite impulse response (FIR) filters Given an input signal X(z) and an output signal Y(z) its transfer function is defined by m m zbzbb zX zY zH −− +++== ....... )( )( )( 1 10 (3) where z stands for frequency in the discrete domain, in difference from Ω , which was frequency in the continuous domain. These filters depend only on the previous and present values of the sampled input function, so the output is typically a weighted average of the values of the input signal. For this reason, they are also called moving average filters. Such a tranfer function only contains zeros, not poles (singularities), so it becomes a very stable approach which, in addition, presents a phase change which is locally linear, so their group delay is assumed as constant. The linear phase change is a consequence of the typical design procedure of these filters, where the Fourier transform of the selected window of application is convoluted to the impulse response expected for the filter. Their main drawback is the finite response they offer to an impulse response, which results in a limited frequency band in the output signal. b) Infinite impulse response (IIR) filters This family involv es tw o dif f erent types o f filters. The first ty pe are the so-called auto-regressive (A R) filters, described by a tranfer function such as Please verify that (1) all pages are present, (2) all figures are acceptable, (3) all fonts and special characters are correct, and (4) all text and figures fit within the margin lines shown on this review document. Return to your MySPIE ToDo list and approve or disapprove this submission. 5856-83 V. 1 (p.3 of 11) / Color: No / Format: A4 / Date: 2005-03-29 13:19:04 SPIE USE: ____ DB Check, ____ Prod Check, Notes: n n zazazX zY zH −− +++ == .....1 1 )( )( )( 1 1 (4) so their transfer function only present a set of poles, without any zero. The stability condition for these filters requires the poles of the function to be inside the unit circle of the root location diagram 7 . This condition needs be v erified along the design. This kind of filters are recursive, as their output depends both on the input value under consideration and on previously calculated output values. The second type of IIR filters manage parts of the FIR and IIR approaches already presented, with transfer functions described by n n m m zaza zbzbb zX zY zH −− −− +++ +++ == ......1 ...... )( )( )( 1 1 1 1 0 (5) T his is t he family usually implemen ted, as far as i t is pa rt ly a moving average filter , a nd pa r tly an auto regressive filter , so they are called ARMA (autoregressive moving average) filters. Their stability needs to be evaluated along the design procedure, as now both zeros and poles are present in the transfer function, but their response to an impulse is not frequency-limited as in the case of FIR filters. This type of filters may be expressed, through a mathematical transform called the z-transform 5 , into an equation in differences such as )(....)1()(...)1()()( 11 0 nkyakyamkxbkxbkxbky n m −−−−−−++−+= (6) where k denotes the order number of the sample in the spatial domain. In the fringe processing approach we are interested in, a FIR approach might suffice, although it would be nothing but the classical smoothing procedure applied to each linear signal. However, whenever a FIR filter is aplied to a given filter frequency specification (defined such as Fig.1), the equivalent IIR filter will always have a lower order, which directly relates to a reduction in the required computing time. The main drawback of these filters is its lack of linerarity in the introduced phase change, which does not allow a prediction o f the position of the fringes and will require an alternative strategy for recovering the original phase after filtering. 3. DESIGN OF THE MATCHED FILTER 3.1 Selection of parameters As presented in fig.1, the main parameters involved in the definition of the filter are its order N and the ones described in fig.1, with frequencies ω p and ω a being equivalent to Ω p and Ω a in the continuous frequenc y d omain; δ p and δ a keep the same meaning of fig.1. The order of the filter becomes determined by the remaining four factors through the expression p a a p N ω ω δ δ log 110 110 log 2 /1 1.0 1.0       − − = − − (7) In our case, available images have individual rows with an extension of 1200pixels. This means the Nyquist frequency of the images is 600 fringes. As far as a low-pass filter is being designed, we have estimated the maximum fringe density in our fringe patterns for different samples of application, which was fixed at 55 fringes. This frequency is to be normalised Please verify that (1) all pages are present, (2) all figures are acceptable, (3) all fonts and special characters are correct, and (4) all text and figures fit within the margin lines shown on this review document. Return to your MySPIE ToDo list and approve or disapprove this submission. 5856-83 V. 1 (p.4 of 11) / Color: No / Format: A4 / Date: 2005-03-29 13:19:04 SPIE USE: ____ DB Check, ____ Prod Check, Notes: to the Nyquist frequency of our fringe-patterns, yielding a normalised frequency of 0.091, which is fixed as the cut-off frequency of our filter . Notice, however, t hat we are speakin g of norma lised fri nge densi ty, so if the fri nges appea r i n one region of the CCD the corresponding cut-off frequency would need to be compensated. Imagining n fringes are the maximum number of fringes which appear, but that they appear along only N used pixels , the corresponding cut-off frequency needs be corrected through Total pixelsfringes p N yquist used pixels N n wnN = (8) yielding a cut-off frequency of 0.071 for n fringes =25 and N used pixels =700pixels, considering the abovementioned n Nyquist value of 600 fringes. This operation could be looked at simply as a correction on the Nyquist frequency depending on the number of pixels registering visible fringes. The remaining parameters are fixed empirically trying to balance the requirement of a very good attenuation for frequencies over ω a , against the requirement of a low-order filter allowing a fast processing of the fringe pattern. Fortunately, fringe-pattern processing does not need abrupt transitions between the attenuated and passing bands. The final filter was designed with a ω a =0.35, a δ p =-3dB, and a δ a =-25dB, yielding a 1:18 reduction in the signal of the highfrequency components. This yields a value of N=2.12, showing a simple filter of order 2 suffices for our purposes. 3.2. Transfer function in the continuous domain Once the parameters have been fixed, the continuous domain transfer function may be designed. One of the most used approaches for this step is the so-called bilinear transform 6,7 , which stablishes a bi-directional correspondence between the continuous and discrete frequency spaces through 21 1 sample z sTz − =+ (9) where T sample is the spatial period of the sampling of the signal in our case. A direct consequence o f using this transform is a compression in the frequency domain, which changes from an infinite continuous frequency axis s to a limited discrete frequency axis z which ranges from 0 to π . This changes the normalized cut-of f frequency of the filter from 0.09 to 0.35 while the attenuation values δ a and δ p are kept constant. This yields the transfer function in the continuous domain defined by 1225.0495.0 1225.0 )( )( )( 2 ++ == ss sX sY sH (10) which corresponds to a low-pass filter under the Butterworth approximation. From eq.10 it becomes evident that the roots of the polinomial (which will become the poles of the filter) have a negative real part, so the stability of the filter response in the continuous domain bec omes ensured. This aspect will be checked again when the response is taken back to the discrete frequency domain. 3.3. Discrete domain implementation Once the filter is designed and validated in the continuous domain, its discrete domain equivalent is obtained by a new application of the bilinear transform (9). Due to the mathematical properties of a double bilinear transformation, the filter dependence on T sample cancels. The result, which finally corresponds to the filter to be implemented, is 2 1 21 628.0517.11 024.0048.0024.0 )( )( )( −− −− +− ++ == zz zz zX zY zH (11) Please verify that (1) all pages are present, (2) all figures are acceptable, (3) all fonts and special characters are correct, and (4) all text and figures fit within the margin lines shown on this review document. Return to your MySPIE ToDo list and approve or disapprove this submission. 5856-83 V. 1 (p.5 of 11) / Color: No / Format: A4 / Date: 2005-03-29 13:19:04 SPIE USE: ____ DB Check, ____ Prod Check, Notes: To properly evaluate the final response of the filter, its response regarding modulus and phase is evaluated. Fig.2 presents such responses, where the upper modulus plot shows the desired attenuations of -3dB at ω p =0.09 and -25dB at ω a =0.35. Notice the slope of the filter is not sharp, due to the low-order used, although it will be shown to succesfully Figure 2: Modulus (up) and phase (down) of the discrete domain filter in eq.11. manage fringe patterns. The lower phase plot corresponds to a second order low-pass filter, showing the non-linearity o f the phase shift which presents an inflexion point near the cut-of f frequency, with a total phase change around 180º. This phase change depending on the frequency will introduce a group delay dependent on the normalized frequency. Figure 3 shows the group delay (left) and the position of the poles of the filter (right) in the discrete domain for the filter described in eq.11. The group delay plot shows the non-linearity o f the filter phase shown in fig.2 introduces a delay o f 4 to 5 samples in the phase of the signal. For a fringe pattern analysis, this means the fringe becomes displaced 4 to 5 pixels from its original position. This drawback will be overcome through filtering each line of the signal in opposite directions: the displacement in the “right to left” filtering becomes compensated with the “left to right” one. The effect of the phase delay becomes then cancelled. This approach is quite common in signal processing, and is usually implemented in general signal processing software 9 .The position of the poles of the filter, in the lower plot of fig.3, shows how both poles of the transfer function in eq.11 (marked with crosses) are placed within the unit circle marked with the dotted line, which is the condition stablished for filter stability in signal processing theory. Finally, the discrete domain transfer function of the filter should be converted to an equation in differences through the ztransform mentioned in eq.6., yielding the final equation of the filter to be implemented as )2(613.0)1(517.1)2(024.0)1(048.0)(024.0)( −−− + − + − + = kykykxkxkxky (11) Notice c oe f ficients o f the transfer function stay the same o f eq.(11), while the exponents of the z term have b een turned into the ordering of the spatial samples of the input X(z) function, turning z -1 into x(k-1), the sample before the one being Please verify that (1) all pages are present, (2) all figures are acceptable, (3) all fonts and special characters are correct, and (4) all text and figures fit within the margin lines shown on this review document. Return to your MySPIE ToDo list and approve or disapprove this submission. 5856-83 V. 1 (p.6 of 11) / Color: No / Format: A4 / Date: 2005-03-29 13:19:04 SPIE USE: ____ DB Check, ____ Prod Check, Notes: considered in the input function. In a similar manner, y(k-2) stands for a sample of the output function of the filter taken two samples before the one being considered, so implying recursivity of the output. Figure 3: Group delay (right) and position of poles (left) of the discrete domain filter in eq.11. 4. APPLICATION TO COMPLEX FRINGE-PATTERN PROCESSING 4.1 Details on implementation The equation in differences is easily implemented into software applications of any given platform, taking into account the original signal (a line in the fringe pattern orthogonal to the fringes) will usually be described as an integer or character variable. x(k) is converted into a floating point value and operated with y(k), which necessarily must be kept as a floating point value to avoid undesirable truncation errors. As far as the final image needs be expressed in integers, a final conversion to 8-bit integers is performed. Finally, once the filtering is performed, we have selected an application where the detection of the maxima of the fringes is performed, although phase extraction procedures could also be equally applied. This has been preferred as it is a very simple and robust approach, although more sophisticated approaches are described elsewhere 2 . We simply detect the point where the local slope of the filtered signal changes its sign, and consider the maximum of the fringe the sample with a higher value. To give further robustness to the system, that point will only be considered a maximum if in the following data the signal decreases a 20% from the value at the maximum. Figure 4: Detection of maxima in filtered signal. Between samples 2 to 4 a change in slope is detected and 3 is considered the maximum as it has the larger signal level. Notice the signal has been already filtered. Please verify that (1) all pages are present, (2) all figures are acceptable, (3) all fonts and special characters are correct, and (4) all text and figures fit within the margin lines shown on this review document. Return to your MySPIE ToDo list and approve or disapprove this submission. 5856-83 V. 1 (p.7 of 11) / Color: No / Format: A4 / Date: 2005-03-29 13:19:04 SPIE USE: ____ DB Check, ____ Prod Check, Notes: 4.2 Application to complex fringe-patterns The filter will now be applied to a complex fringe-pattern obtained through a digital Ronchi deflectometry set-up, where the wavefront crosses a progressive addition ophthalmic lens before passing through a low frequency test with square transmittance, yielding the fringe pattern depicted in fig.5. The size of the image is 1200x1170 pixels, and it is acquired from behind a diffusiv e screen, so an important noise level is present in the image. Notice also the different illumination levels in different regions of the image, both in the fringes and in the bias, which make a binarisation of the image impossible without very important losses of information. The distance vision area of the progressive lens is placed on the upper side of the image, while the near-vision area is on the lower side. The power progression is quite evident. As far as the total image is composed after different acquisitions, a noticeable change in average intensity may be appreciated in the middle of the image, adding further complexity to the processing. It can also be appreciated how the signal is saturated in the sides of the progression area. Figure 5: Complex wavefront to be processed As explained, the signal in each row will be processed separately. One of the key aspects which needs be confirmed is the negligible group delay introduced in the signal when filtering twice along opposite directions is applied. This has been analysed in three different regions in image, corresponding to the distance vision area (row 250), the power progression area (row 650) and the near vision area (row 1100), and is presented in figure 6. The selection of rows corresponds to quite different signals allowing to show the ro bustness o f the filtering strategy propos ed, which performs without any modification in the three cases. For each row (labelled as a to c), figure 6 presents the non-processed data (labelled as 1), the data after the first filtering (labelled as 2), and the final output data after two passes of the filter in opposite directions (labelled as 3). An enlarged plot with these three signals superimposed (continuous line for the original signal, large dash for the single-filtered signal and short dash for the double-filtered signal) for a typical fringe is also included and labelled as 4 in all cases. It may be seen in each plot how the group delay introduced is negligible after filtering twice in opposite directions, regardless the complexity of the signal. However, notice the group delay of some pixels visible in the single filtered signal (long dash) which disappears in the double-filtered signal (short dash). The variety of situations presented in these three plots includes different number of fringes (from 28 fringes in fig.6b to 15 in fig 6c), oscillations in the maximum and minimum signal level in the fringes within the signal (all cases), local bias levels of the signal bigger than the signal level at some maxima (fig.6b), saturated signals (fig.6b), and quite regular, smooth patterns (fig.6c). The expected highfrequency reduction in the signal is also evident in all cases. The designed filter perform equally well in all situations, proving its robustness under different signal conditions. Please verify that (1) all pages are present, (2) all figures are acceptable, (3) all fonts and special characters are correct, and (4) all text and figures fit within the margin lines shown on this review document. Return to your MySPIE ToDo list and approve or disapprove this submission. 5856-83 V. 1 (p.8 of 11) / Color: No / Format: A4 / Date: 2005-03-29 13:19:04 SPIE USE: ____ DB Check, ____ Prod Check, Notes: (a) (b) (c) Figure 6: Original and processed patterns for a) Row 250, up; b) Row 650,middle; c) Row 1100, below. See text. Please verify that (1) all pages are present, (2) all figures are acceptable, (3) all fonts and special characters are correct, and (4) all text and figures fit within the margin lines shown on this review document. Return to your MySPIE ToDo list and approve or disapprove this submission. 5856-83 V. 1 (p.9 of 11) / Color: No / Format: A4 / Date: 2005-03-29 13:19:04 SPIE USE: ____ DB Check, ____ Prod Check, Notes: