Surface characterization by structure function analysis

The structure function is a tool for characterizing technical surfaces which exhibits a number of advantages over Fourier-based analysis methods. So it is optimally suited for analyzing the height distributions of surfaces measured by full-ﬁeld non-contacting methods. After the deﬁnition of line-and area-structure function and offering effective procedures for their calculation this tutorial paper presents examples using simulated and measured data of machined surfaces as well as optical components. Comparisons with the results of Fourier-based evaluations clearly prove the advantages of structure function analysis.


INTRODUCTION
The last decades have seen a rapid development of optical noncontacting metrology systems for measuring surface contours of technical components with high precision in the microscopic as well as in the macroscopic realms [1,2].For microscopic components the confocal scanning methods are state of the art [3].For macroscopic objects most methods are based on the triangulation principle such as photogrammetric methods [4].If patterns are projected, these may be points, as in the Shack-Hartmann methods [5], lines, grids, or even complex synthesized patterns in reverse engineering applications [6].Fringes are projected onto the surface in Fouriertransform profilometry [7] or moire-techniques [8].Specularly reflecting surfaces are measured by deflectometry [9,10].Interferometric methods like holographic contouring also find numerous application fields [11].The resolution achievable by these methods can be further improved by using multiple redundant images employing phase-shifting [12] or Graycode [13].
A common problem of all these metrologic methods is the handling of the raw data [5].Normally the measured surface data are delivered as a point cloud which is stored in a suitable format in computer.But the utilization of the data for quality assurance purposes demands a representation of the data in a form appropriate for the problem to be solved.The questions to be answered by evaluation of the raw data can be manifold: In quality control one may be interested in the existence of defects such as holes, dents, bumps, scratches, or other flaws which normally manifest as localized high spatial frequency variations of the surface heights.On the other hand there may be interest in global criteria like e. g. periodicities, waviness, lay, or roughness.
Often Fourier-transform based methods are the first choice [14,15].So a few Fourier descriptors can indicate the existence of defects, or the power spectral density can describe global parameters of the measured surface data.Nevertheless the Fourier-approach exhibits severe drawbacks: Discontinuities at the edges of the sampling interval in 1D or 2D spread power all across the spectrum, an effect which is minimized by the application of window functions, like the well-known Hanning window.But this has the undesirable consequence that the measured samples enter the calculation of the spectra with unequal weights.Furthermore the 2D Fourier-transform assumes rectangular fields.In practice we often have arbitrary aperture shapes, so we can use zero-padding or choose a sub-aperture and assume it being sufficiently representative of the true surface characteristics.As a result we insert extrapolated data not really measured or discard measured values.
If the intention of our measurements is defect detection, a wavelet approach [16] was suggested as an alternative.If on the other hand one is more interested in the global criteria, the structure function is a suitable tool and will be presented in more detail in the following.

DEFINITION OF STRUCTURE FUNCTION
The structure function originally was introduced by Kolmogorov in 1941 for analysis of statistical problems associated with turbulence theory [17,18].In the meantime it found numerous applications in diverse fields: It was used for the characterization of technical surfaces [19], of optical surfaces [20,21], of velocity fields in turbulent flows [17,18,22], and of phase and frequency instabilities in frequency sources [23], as well as for time series analysis in astronomy [24,25], to name just a few.
Let f ( x) be the measured data, which may be the height profile of a rough surface [19] or an optical surface [20,21], a flow velocity [22], optical flux [25], optical phase differences [26] or similar data sets.Generally f ( x) can be treated as a random process.Then the structure function of order d is defined by [22].Most important and the only one used in the following is the structure function of order 2, so the subscript d will be omitted: S = S d .
In detail, we have in one dimension the so called line structure function If the x = (x, y) is two-dimensional, as e. g. the coordinates of a surface with heights f (x, y), then the 2D structure function, also called area structure function, is To allow a meaningful comparison of the structure function in various spatial frequency bands, the normalized structure function NSF(x) can be introduced [25] NSF{ f }(x ) = S(x ) There are applications where the 1D line structure function suffices for evaluation of 2D data sets.Then each pixel pair is not represented by the pair of distances in each coordinate direction (x 2 − x 1 , y 2 − y 1 ) but only by the distance (x 2 − x 1 ) 2 + (y 2 − y 1 ) 2 .This already allows separation into characteristics like figure, roughness, and mid-spatial frequencies [21].
As an example we analyse a Zernike polynomial Z

CALCULATION OF THE STRUCTURE FUNCTION
In practice we are dealing with discrete measured data { f i , i = 1, . . ., I}.Then the structure function is calculated as with N < I. Due to the squaring we have an even structure function S(−n) = S(n) for all n and also S(0) = 0, so it is not necessary to calculate redundant information and it suffices to determine the structure function only for positive differences n In two dimensions with { f i,j , i = 1, . . ., I, j = 1, . . ., J} we calculate the finite sum with N < I and M < J.This area structure function exhibits point symmetry S(−n, m) = S(n, −m) for all n and m and also S(0, 0) = 0 so it suffices to calculate the area structure function in two quadrants.The parameters N and M are in the same range as the I and J and theoretically can be chosen up to N = I − 1, M = J − 1.But then for large n or m we average over only a few squared differences, so the determined values become unreliable due to noise and roughness.Therefore a choice of N and M as at most 80 % to 90 % of the magnitude of I and J is recommended.
It has to be mentioned that the calculation of S(n, m) is not separable into a product of two sums, therefore the computational effort in the two-dimensional case is proportional to I 2 J 2 multiplications.In the case of 2D data sets to be analyzed with data value numbers in the range I × J = 10 3 × 10 3 this would lead to a number of floating-point multiplications proportional to 10 12 /4.A significantly lower number of multiplications is required by the following algorithm that also works for arbitrarily shaped apertures.Let I = max{i} and J = max{j} be the maximum numbers of points in iand jdirection and assume a rectangle of I × J measured values, that contains all actually measured points.Then the algorithm works as follows: • The rectangle of the I × J measured pixels is divided into I × J subareas, each of the same of length about 10 × 10 measured values.
• For each of these subareas we initialize its sum by zero and also its count by zero.
• Using a random number generator we select a random pair of points.If one of these points falls outside the aperture, the random selection is repeated.
• The distance vector r between these two points is calculated.
• We calculate the squared difference of the measured values at the two points.
• The subarea containing vector r gets its sum increased by the just calculated squared difference and its count increased by one.
• The steps beginning with the random choice of a point are repeated about 10 5 to 10 6 times.
• Finally for each subarea its cumulated sum is divided by its cumulated count.
This algorithm determines an area structure function which approximates the averaged exact area structure function, where averaging is performed over rectangles of size (I/I ) × (J/J ).Due to the random choice of point pairs, small differences will occur more frequently.Thus the approximation is highly reliable for differences close to zero, but becomes noisy when approaching large differences, requiring both randomly chosen points to be positioned at opposite ends in the measurement aperture, see results in Section 5.
The big advantage is the significant decrease of the number of multiplications which is of the order of 10 3 to 10 4 , see Section 5.

PROPERTIES OF THE STRUCTURE FUNCTION
For statistically stationary processes the structure function S( x) of order 2 is related to the autocovariance function . It contains the same information as the autocovariance function and its Fourier transform, the power spectral density function, but offers some practical advantages: it is stable and easy to compute, it does not impose a periodogram model on the measured surface, and it avoids the singularity at the origin of the autocovariance function [19].Also the structure function can be related to the autocorrelation function φ f f (x ): Let σ 2 be the variance of the underlying stochastic process then see Appendix A. Although the structure function thus carries principally the same information as the power spectral density function, it does not have the disadvantages of the Fourier-transform based indicators named in the Introduction.The structure function can easily deal with arbitrary apertures and does not need any windowing.A more general approach to the structure function is given in Appendix B.
The structure function has a minimal value of 0 at pixel distance zero, S(0) = 0.This corresponds to the maximum of the autocorrelation function at zero, φ f f (0) = 0.The squaring in the calculation of S 2 causes the evenness of S 2 , so the 1D structure function S 2 (x ) normally is only used for positive pixel distances x ≥ 0.
The advantages of the structure function approach compared to Fourier-based methods are pointed out in the next example, Figure 2. We assume an irregularly shaped surface which furthermore shows a hole, Figure 2(a).In practice it does not matter whether this hole physically exists or only represents corrupted data which are to be excluded from evaluation.In Figure 2(a) bright area represents surface data points, dark stands for no data.This results in possible and impossible point pair distances.and meaningful approximation, at least at all point pairs not too close to the margins.Figure 7 clearly demonstrates the decreasing noise when increasing the number of randomly chosen pairs.The long spatial wavelength components in vertical direction in the measured surface can be recognized in the area structure function of Figure 7.The appearance of parallel fringes in the area structure function gives strong evidence that the detected waviness is a global one, it does not vary from one column to another.There is no direct correlation between the local amplitudes of the area structure function and the height of the waviness in the measured data.But due to its integrating nature the structure function will detect regular patterns like lay or waviness even if these are embedded in noise, which may be generated by surface roughness or speckles, if coherent light was used for measurement.Furthermore in the interpretation of Figure 7 one has to pay attention to the subarea size: the numbering along the axes relates to this subarea size.In our example we have chosen subareas of 10 × 10 pixels, so the large differences at distances 100 and 275 indicated in Figure 6 are reflected at subarea distances of 10 and 27.5.

APPLICATION OF THE STRUCTURE FUNCTION
The next example is from characterization of optical components.Here usually the transmitted optical wave fronts are investigated, so we must analyse the structure of optical phase distributions.If the exit pupil of an optical system to be tested is circular, which is the case in most applications, then the aberrations present in the system can be represented in terms of Zernike polynomials.These are orthogonal and normalized within a circle of unit radius [27].Then the phase error φ A (x, y) is represented by an expansion into the Zernike polynomials Z k (ρ, φ) where ρ is the radial coordinate within the unit circle and φ is the polar angle.Each Zernike polynomial is of the form and the root-mean square error φ A due to aberrations is As He et al [21] pointed out, the residual error in φ A (ρ, φ) can be determined by subtracting the first K Zernike terms, then this error should be analyzed by the area structure function for detecting e. g. ripples.7 (ρ, φ) is odd in the vertical direction, its corresponding area structure function by principle is even in the horizontal as well as in the vertical direction.He et al [20,21] have pointed out the capability of the area structure function for detecting surface anisotropy and distinguishing different error geometries when investigating measured phase distributions, while the line structure function does not exhibit this ability.
Figure 9 gives the last example, the height distribution along an airplane wing section measured by fringe projection.We observe that the wing section does not fill the full rectangular area of 919 × 481 measurement points.Especially there is a hole which is clearly seen in the mask image of Figure 10 with all measurement points marked white and the points with no measured value marked black.Here we have a typical example of a non-rectangular aperture.First we calculate all 1D line structure functions along the columns, see Figure 11.The columns are numbered from 1 to 919 along the abscissa, while along the ordinate we have spatial differences, the arguments of the structure functions.For the low numbered columns the structure function values of large differences are lacking, because due to the non-rectangular aperture there do not exist Only positive spatial differences.
two measured points far apart.The same effect, but not as pronounced as before, is seen for the high numbered columns.
Local minima appear at differences around 130, 260 and 390 pixels, which indicates a waviness with a periodicity around 130 pixels.This also is made obvious in the average of all these line structure functions displayed at the left of Figure 11.The outliers at point differences higher than 420 are due to the fact that there is averaging over only very few point pair differences.
The 1D line structure functions along all lines are shown in Figure 12.Here no distinct local minima are present in the structure functions, neither in the integral over all line structure functions displayed above, nor locally in the individual line structure function.Thus no waviness in horizontal direction is detected.Nevertheless we observe increasing amplitudes with increasing differences (from right to left) in the integrated structure function.This indicates a slight overall tilt of the measured structure, which can be eliminated before the structure function calculation, see Section 6.The area structure function calculated using 5 × 10 6 randomly chosen point pairs is given in Figure .13. Especially the average structure function gained by integration over all lines or columns, confirm the aforementioned statements regarding the waviness.The area structure function in Figure 13 is approximately point symmetric.It is not perfectly point symmetric, since due to the stochastic choice of point pairs not for each point pair ((x 1 , y 1 ), (x 2 , y 2 )) also the corresponding symmetric pair ((x 2 , y 2 ), (x 1 , y 1 )) is used for calculation.This procedure maps each point difference into the corresponding one of the four quadrants.A better procedure would be the calculation of the 2D structure function in only two quadrants, thus using in the average twice as many samples at each point difference of the structure function.

STRUCTURE FUNCTION AND FILTERING
Generally the form of the structure function is influenced by an uneven background in the data set.In many cases the measured underlying background of a plane surface is a linear function due to a tilt of the surface, a parabolic function indicating a bending, or a bellshaped curve caused by irregular illumination when using optical intensity based measurement methods, to name just a few.The effect of determining such background by linear regression and removing it before calculating the structure function is shown in Figures 14  to 16.The measured heights of a rectangular surface are displayed as gray values in Figure 14.The heights along a central line are given in Figure 15.A cubic function describing a global deformation is fitted by linear regression, Figure 16(a), the differences of the heights from this parabola are shown in Figure 16(b).The structure function calculated from the original values is displayed in Figure 16(c), while the structure function of the heights without background can be seen in Fig. 16(d).The dissimilarity between Figures 16(c) and d can be interpreted as follows: The difference between the first left maximum and the rightmost local minimum is amplified by the uneven background.This difference is reflected in the structure function as a local maximum near a 360 pixel shift, the right local maximum in Figure 16(c).In contrast, the left local maximum is the higher one in Figure 16(d), resembling the differences between neighboring opposite local extrema in Figure 16(b).The measured values exhibit a waviness with a period of about 260 pixels.This is seen in Figure 16(d) as a central local minimum.As a consequence one must conclude that a previous background correction yields more reliable structure functions and thus is strongly recommended.
The background elimination can be viewed as a sort of highpass-filtering.It has to be applied judiciously depending on the spatial frequencies of interest.Furthermore one has to keep in mind that each value of S(n) is an average over I − n differences or a subset thereof.So for rough surfaces the values S(n) for small n are more reliable than those for large n.This effect can be diminished by a low-pass-filter with a cutoff-frequency chosen to eliminate the high-frequency surface height fluctuations due to roughness.

CONCLUSION AND OUTLOOK
In this paper we have defined and analyzed the structure function and have shown its ability to characterize surfaces whose height distribution is measured by one of the numerous optical noncontacting methods.Various ways to extract 1D structure function information from the 2D height distribution have been demonstrated, some of them -those aver-aging pointwise over all line structure functions -for the first time.The advantages of the structure function compared to Fourier-based approaches have been pointed out and shown by example.Especially if the contour of the measured surface is not rectangular or if the continuation at the surface's margins is not smooth, then the structure function approach has to be favoured, since it exhibits no leakage nor aliasing.The structure function's applicability is demonstrated by simulated as well as practically measured examples like surface height variations or optical phase distributions.We have seen that the positions of local extrema indicate precisely the existence or non-existence of waviness as well as its spatial frequency.What is still lacking, is a general procedure to draw significant parameters from the structure function which characterize the investigated surface, such as indicating and quantifying lay, waviness, slope, roughness and other parameters.
Although the emphasis of the paper was on evaluation of optically measured height distributions, the capability of the structure function approach to phase error characterization in the testing of optics also was mentioned.The results regarding the structure function presented in this paper should draw the attention of those analyzing technical or optical surfaces to structure function analysis and its manifold capabilities.
The process g(t) has a stationary Mth increment, if the following averages exist for all real T and τ and do not depend on the instant t: ∆ (M) g(t; τ) • ∆ (M) g(t + T; τ) = D (M) g (T; τ) (14) In other words, the Mth increment has a time independent mean and its autocorrelation depends only on the time difference T (wide sense stationarity) [23].
By definition D So the 1st structure function of the process is D g (τ) = [∆ (1) which coincides with our previous definition of the structure function.

Figure 2 (
b) displays the first quadrant of the plane of point pair distances indicating possible pairs in bright, impossible pairs in dark.E. g. for the distance (400,100), represented in red, no valid point pair can be found in Figure2(a).Three point pairs with distance (400,100) are given in Figure2(a), none with both endpoints can be found inside the surface.On the other hand point pair distance (90,250) shown in green is a possible one, also three point pairs for this distance are drawn into Figure2(a).The simulated measurement values describe a cosine-valued displacement A cos( f πx) with A = 1.5, f = 0.038 and the x-coordinate along the horizontal axis, Figure2(c).The resulting area structure function can be seen in Figure2(d).It reflects the periodic nature of the displacement function in horizontal direction as expected.Only at the margins of the allowed area can higher fluctuations be seen.For these distances only a few pixel pairs can be found in the measured data, so errors are not compensated by averaging a high number of measured points.For comparison purposes a Fourier based evaluation was performed.The pixels outside the surface were given the value zero to fill the full 512 × 512-pixel square -the common zero padding.The Fourier transform F(u, v) = F { f (x, y)} of this real distribution is calculated employing the FFT-procedure and the power spectral density S(u, v) is derived pointwise by S(u, v) = F(u, v)F * (u, v).We then obtained the autocorrelation φ f f (x, y) by applying the inverse Fourier transform φ f f (x, y) = F −1 {S(u, v)} where the x, y can be interpreted as shifts or point pair distances.1.0 − φ f f (x, y) is given in Figure 2(e), it should be proportional to the structure function, Figure 2(d).For better visual comparison this distribution of Figure 2(e) is masked by the values 1. and 0. of Figure 2(b) and in the resulting Figure 2(f) we recognize severe differences to Figure 2(d), which result from aliasing and zero padding in the Fourier approach.The advantages of structure function evaluation are therefore obvious.

Figure 3 FIG. 2 FIG. 3 2 FIG. 4 10 FIG. 5 FIG. 6 FIG. 7
Figure 3 shows an example from practice, a measured height profile of 445 × 403 points.As can be seen in the gray-scale display, only values in a rectangular, inclined aperture are measured.The height-values along a single column, here arbitrarily the one with index 201, are extracted, see Figure 4. Along this line the 1D structure function is calculated with parameters I = 445 and N = 400, see Figure 5.This line structure function starts with zero at the origin and has a local minimum at a displacement of about 175 pixels, and two local

Figure 8 (
Figure 8(a) shows the Zernike polynomial Z 1 7 (ρ, φ) and Figure 8(b) displays its area structure function.While the Zernike polynomial Z 17 (ρ, φ) is odd in the vertical direction, its corresponding area structure function by principle is even in the horizontal as well as in the vertical direction.He et al[20,21] have pointed out the capability of the area structure function for detecting surface anisotropy and distinguishing different error geometries when investigating measured phase distributions, while the line structure function does not exhibit this ability.

FIG. 10 FIG. 12
FIG. 10 Mask image.White: measured point; black: no measurement FIG. 13 2D-structure function and integrated values along all line (left) and column coordinates (above).Abscissa and ordinate: spatial differences.Positive and negative differences of spatial positions in each direction.