Intelligent sampling for the measurement of structured surfaces

Uniform sampling, sometimes referred to as uniform stratified sampling, allocates sampling points in a regular latticed pattern, for example, the most classical regular 'square grid' pattern which is investigated in this study. Investigations on different uniform sampling methods are described elsewhere [31, 32]. The 'square grid' uniform sampling method is implemented by allocating a random start point ${\bm p}_{0} = [x_0 ,y_0 ]_{x_0 ,y_0 \in \mathbb{R}}$ in a plane and subsequent periodic translations $d = [d_x ,d_y ]_{d_x ,d_y \in\mathbb{R} }$ and duplications on the x/y directions. Thus, a 'square grid' pattern P can be constructed (see figure 3(a)):

$$\begin{equation*} {\bf P} = \{ {\bm p}_{i,j}: {\bm p}_{i,j} = {\bm p}_0 + [id_x ,jd_y ],i,j \in \mathbb{Z}\}. \end{equation*}$$

Zoom In Zoom Out Reset image size

Jittered uniform sampling is a simple variation of uniform sampling. Simply added with a random position jittering sequence $\Delta = \{ \Delta _{i,j} \in \big( - \frac{\bm d}{2},\frac{\bm d}{2}),i,j \in\mathbb{Z} \}$, a uniform sampling pattern becomes a jittered uniform sampling pattern (see figure 3(b)):

$$\begin{equation*} {\bf P}_{{\rm Jittered}} = \{{\bm p}'_{i,j} :{\bm p}'_{i,j} = {\bm p}_{i,j} + {\rm \Delta }_{i,j} {\rm , }{\bm p}_{i,j} \in {\bf P},{\rm \Delta }_{i,j} \in {\rm \Delta },i,j \in {\mathbb Z}\} . \end{equation*}$$

By optimizing the sampling pattern, an estimation error of the total or average properties of the population can be reduced. Two classical low-discrepancy patterns—Hammersley and Halton pattern—are introduced here by considering their two-dimensional case:

$$\begin{equation*} \begin{array}{*{20}c} {{\bf P}_{{\rm Hammsley}} = \{ {\bm p}_i :{\bm p}_i = [i/N,\Phi _2 (i)],{\rm }i \in\mathbb{Z} \} } \\ {{\bf P}_{{\rm Halton}} = \{ {\bm p}_i :{\bm p}_i = [\Phi _2 (i),\Phi _3 (i)],{\rm }i \in\mathbb{Z} \} } \\ \end{array}, \end{equation*}$$

in which Φ_b(i) is the radical inverse function [5] Φ_b(i) = 0.d₁d₂...d_m, where d_i satisfies $\sum\limits_{k = 1}^\infty {d_k b^{k - 1} } = i,{\rm }d_i \in \{ 0,1, \ldots ,b\}$. The radical inverse function converts a non-negative integer i into a floating-point number within [0, 1). For example, the first ten numbers of Φ₂ are $\left[ {\frac{{\rm 1}}{{\rm 2}}{\rm ,}\frac{{\rm 1}}{{\rm 4}}{\rm ,}\frac{{\rm 3}}{{\rm 4}}{\rm ,}\frac{{\rm 1}}{{\rm 8}}{\rm ,}\frac{{\rm 5}}{{\rm 8}}{\rm ,}\frac{{\rm 3}}{{\rm 8}}{\rm ,}\frac{{\rm 7}}{{\rm 8}}{\rm ,}\frac{{\rm 1}}{{{\rm 16}}}{\rm ,}\frac{{{\rm 15}}}{{{\rm 16}}}{\rm ,}\frac{{\rm 5}}{{{\rm 16}}}} \right]$. The Hammersley and Halton patterns of 100 points are illustrated in figures 3(c) and (d), respectively.

Three adaptive sampling methods have been studied in this research and these comprise the following: sequential profiling adaptive sampling, triangle patch adaptive subdivision sampling and rectangle patch adaptive subdivision sampling. Adaptive sampling has no fixed sampling pattern or mathematical expressions similar to the four methods above. Adaptive sampling methods are expressed by ordered algorithm procedures.

3.4.1. Sequential profiling adaptive sampling

Sequential profiling adaptive sampling [6] has been developed based on Shih's indirect sampling [10]. Considering the raster scanning mechanism used in conventional stylus profilometers, this method is expected to work on current stylus instruments to enable efficient measurements. Sequential profiling adaptive sampling is a non-CAD model-based method which consists of a core algorithm—profile adaptive compression—and works in two sequential stages.

The profile adaptive compression algorithm requires an initial profile scanning as a reference which employs a high-density sample size setting, e.g. the instrument allowed highest sample density. The key samples are then selected from the initial samples by recursively examining the reconstruction error compared to a pre-defined threshold. Specifically, in the first stage, an approximate measurement (for example, take several profile scans in the y-direction at random or uniformly selected x-positions, see figure 4(b)) is carried out and the y coordinates of the key scanning positions can thus be determined based on the adaptive compression algorithm. In the second stage, a fine adaptive sampling is implemented in the x-direction at each of the key scanning position. With the aid of the adaptive compression algorithm, the key samples of each x-direction scan are collected which constitute the final sampling result. As illustrated in figure 4(c), a final sampling result can thus be achieved. This method has been demonstrated with apparent advantages over the four sampling patterns for the measurement of structured surfaces. A detailed description of this technique can be found elsewhere [6].

Zoom In Zoom Out Reset image size

3.4.2. Triangle patch and rectangle patch adaptive subdivision sampling

Triangle patch and rectangle patch adaptive subdivision samplings are two CAD model-based methods. They require a CAD model or a preliminary measurement (with a simple sample design, e.g. uniform sampling) for the determination of the adaptive samples. The two methods have been introduced as 'direct sampling' [10]; hereby, a brief introduction is given and a minor modification to the error evaluation criteria is made. The algorithm is as follows.

• 1.  
  Triangle patch adaptive subdivision sampling.

  • (1)  
    Select a rectangle region on a surface as the sampling object.
  • (2)  
    Select four initial points on the extreme corners of the rectangle region and group the four corner points into two triangles.
  • (3)  
    Subdivide each triangle into four triangles by inserting three points on the centre of the edges of the triangle, as in figure 5(a).
  • (4)  
    Evaluate the reconstruction error³ of each triangle.
  • (5)  
    If the error is greater than a preset threshold, then repeat steps 3 and 4. Otherwise, stop.
• 2.  
  Rectangle patch adaptive subdivision sampling.

  • (1)  
    Select a rectangle region on a surface as the sampling object.
  • (2)  
    Select four initial points on the extreme corners of the rectangle region.
  • (3)  
    Subdivide the projected area into four rectangles by inserting five points as shown in figure 5(b).
  • (4)  
    Evaluate the reconstruction error¹ of each rectangle.
  • (5)  
    If the error is greater than a preset threshold, then repeat steps 3 and 4. Otherwise, stop.

Zoom In Zoom Out Reset image size

Representative sample patterns generated by the three methods are illustrated in figure 6 in which the ideal surface tested in figure 4(a) is sampled respectively using each method above with a sample size of approximately 1500 points. It is found that the sequential profiling adaptive sampling pattern has no regular distribution of sample points on the feature edges. Triangle patch and rectangle patch adaptive division samplings generate more regularly designed sample patterns. Specifically, the former has dense samples regularly on the feature edges, while the latter yield more dense samples at the feature lateral corners.

Zoom In Zoom Out Reset image size

The tensor product method has been widely used for the reconstruction of regular lattice data (for example, the uniform sampling result in figure 3(a) or partially regularly latticed data in figure 6(a)); this is due to this method's high numerical stability and computational efficiency. The tensor product method presents a surface as a tensor product of two bases, for example,

$$\begin{equation*} \left\{ \begin{array}{@{}c@{}} {\phi (x,{\bm a}) = \sum\limits_{k = 1}^S {a_k \phi _k (x)} } \\ {\varphi (y,{\bm b}) = \sum\limits_{l = 1}^T {b_l \varphi _l (y)} } \\ \end{array} \right.,\end{equation*}$$

in the x and y directions independently. Thus, the surface can be expressed as

$$\begin{eqnarray*} &&\fl z = \phi (x,{\bm a}) \times \varphi (y,{\bm b}) = \sum\limits_{k = 1}^S {\sum\limits_{l = 1}^T {a_k b_l \phi _k (x)\varphi _l (y)} } \nonumber\\ &&= \sum\limits_{k = 1}^S {\sum\limits_{l = 1}^T {c_{k,l} \psi _{k,l} (x,y)} } ,\end{eqnarray*}$$

where {ϕ_k} and {φ_l} are preset base functions and the coefficient vectors a and b should be calculated from the data. Chebyshev polynomials, polynomial splines and B-splines provide base functions for the tensor product surface reconstruction [2]. Considering the smoothness and computational stabilities, second- (linear) and fourth-order (cubic) B-spline [33] basis functions are adopted in this research.

Intelligent sampling always results in non-regular latticed data which cannot be solved by the tensor product method. Triangulation-based methods are a simple and stable substitution. For example, Delaunay triangulation [34, 35] establishes neighbourhood connections among the data points with the Delaunay triangulation algorithm, which neglects all the non-neighbouring points in the Voronoi diagram of the given points and avoids poorly shaped triangles. Following this structuring process, regional reconstructions [35] (linear or cubic) within each triangle patch can be carried out. These methods are able to guarantee a reconstruction to arbitrary accuracy if the sample points are dense enough, which provides the theoretical foundation for developing new reconstruction techniques.

Considering that the amount of sample points for surface measurement is usually large, RBF-based interpolations or fits are not generally recommended. For example, it has been stated elsewhere [36] that RBF-based reconstructions may be very unstable, computationally complex and memory consuming; the RBF-based reconstructions are only employed when data points are no more than several thousand in number.

Typical examples using the tensor product reconstruction methods and the Delaunay triangulation method are presented in figure 7 in which the sample result shown in figure 6(a) is tested. The performance differences of the reconstruction results are clearly shown. Selecting an appropriate method for reconstruction usually depends on many conditions, such as the surface complexity, distribution of the sample points, accuracy and efficiency requirements, and so on. In this study, all the potential methods are tested and the best one (with the minimum residual error) is selected for use.

Zoom In Zoom Out Reset image size

Traditional surface geometry descriptors and the traditional statistical parameters such as Sa, Sq or Ssk [37] have been reported to have a lack of efficiency when used to characterize structured surfaces [4]. In a new feature-based characterization system [17, 38], surfaces are treated as a composition of diverse geometric features. Characterization of each individually recognized feature and a statistic of the concerned feature attributes have been developed. Therefore, the purely statistical methods have not been used in this study although they have been used in earlier research work [8, 9, 14].

Estimations of the statistical values (for example, the RMS value, maximum, etc) of the height residuals are usually employed when measuring freeform surfaces. In this way, the CAD model and the measured surface are compared based on height information, and a residual map can then be obtained. Evaluation of the residual map would seem appropriate forperformance comparison of different sampling methods. This solution has been used in nearly all of the current research in metrology [6, 10–12, 15, 26].

In this study, the RMS height residuals error is employed using the following:

$$\begin{equation*} e = \sqrt {\sum\limits_{i = 1}^N {(z_i - z'_i )^2 } } ,\end{equation*}$$

where z_i are each value of the standard high-density sampled data matrix, z'_i is each value of the reconstructed surface data matrix which has the same sample size as the standard data and N is the size of the high-density sampled data matrix.

Feature-based characterization has been recognized as being of high importance in advanced metrology techniques [39], particularly for the characterization of structured surfaces by extracting micro-scale-dimensional parameters [40, 41]. By sequentially employing the 'areal feature segmentation' [16, 42], 'boundary segmentation' [41, 43], 'dimensional parameter selection' (such as the defined attributes in [17]) and the 'parameter calculation and statistics' [38], the concerned micro-scale-dimensional parameters can be extracted. These parameters are then used for the performance comparison of different sampling methods. For example, the mean absolute deviations of these evaluation parameters from that of the original high-density sampled results are investigated in this paper:

$$\begin{equation*} e = {\rm mean}(\left| {p_0 - p_i } \right|), \end{equation*}$$

in which p₀ is an evaluation parameter (for example, step height, roundness) extracted from the original high-density sampled surface or CAD models, and p_i are the corresponding parameters extracted from a reconstructed surface that was sampled at a lower density. For example, figure 8(a) shows a Fresnel lens surface (table 1, (c)) reconstructed from a 2500 Hammersley pattern sample and the areal segmentation is carried out by extracting the feature edges (bold yellow curves). Evaluation parameters p_i such as radius and roundness of a circle feature can be consequently estimated (figure 8(b)). If repeated tests are carried out, a mean value of the results from the same sampling setting is taken.

Zoom In Zoom Out Reset image size

Table 1. The three typical structured surface specimens and the experimental settings.

a

^aThe tested sample sizes are selected based on the following criteria: (1) the tested samples sizes are representatively selected which indicates they might be normally used in practical measurements; (2) the tested sample cannot be too small in which case severe reconstruction distortion occurs; (3) the tested sample size cannot be too large in which case reconstruction error has minor fluctuations and the evaluation process may be time consuming.

General structured surfaces contain three basic types of feature groupings. These comprise linear patterns, tessellations and rotationally symmetric patterns [4]. Three representative high-precision structured surface specimens were used in this research to validate the performance of different sampling methods. High-precision measurements of the three surfaces are presented in table 1, (a)–(c), which include a five parallel-groove calibration artefact, a nine pits-crossed-grating calibration artefact and a Fresnel lens central patch. These high-density sampled results are used as the references for comparison. A performance evaluation procedure is described as follows.

• (1)  
  Standard data preparation.For a given real surface, obtain the standard measurement result (e.g. a 1024 × 1024 regular lattice data) using the instrument allowed high-density sample setting. The standard-feature-related parameters (e.g. groove-width, step height, see table 1, which are selected in consideration of the main functions of the structured surfaces) are characterized for later uses.
• (2)  
  Sampling.Re-sample the standard surface data using different sampling methods and sample sizes. In this study, seven sampling methods and six different sample sizes (see table 1) on each of the specimens are tested.
• (3)  
  Reconstruction.Potential reconstruction methods are then employed to reconstruct the 'continuous' surface with the same sample design as the originals—the standard measurement result. The best reconstruction results with the lowest RMS height residuals are selected for use.
• (4)  
  Performance evaluation.Extract the RMS height residuals of each reconstructed data from the standard data; extract feature related parameters from each reconstructed surface and the differences from the standard measurement results are calculated. The smaller the height residual or a parameter difference, the better the performance of the sampling design used.

The results in table 2 show the sampling errors—the discrepancies of each evaluation parameter between each sampling-reconstruction result and the standard high-density sample result. It indicates that the sampling error has a power function-like relationship against sample size

$$\begin{equation*} y = cx^b ,\end{equation*}$$

where y is sampling error, x denotes the sample sizes and c and b need to be calculated to give the best fit function. Other researchers have shown similar results [6, 8, 10, 23]. A linear plot of the sampling errors against the sample sizes cannot render a clear performance comparison when sample size increases.

Table 2. Deviations of the evaluation parameters from the standard result for samples 1, 2 and 3 (log–log plots).

Since a linear relationship existed between ln y and ln x deduced from the above power function

$$\begin{equation*} \ln y = \ln c + b\ln x,\end{equation*}$$

log–log plotting is employed in this research thus the sampling performance can be shown evenly. By plotting the sampling-reconstruction errors against the sample sizes and giving the best fitting power functions, the sampling errors of each evaluation parameters extracted from the numerical experimental results are presented in table 2.

The following conclusions can be drawn in sequence.

• (1)  
  Adaptive sampling methods usually have prominent advantages over other methods in terms of minimizing the sampling error (height residuals and feature parameters) for structured surfaces.
• (2)  
  Uniform sampling, jittered uniform sampling, Hammersley pattern and Halton pattern samplings have close capabilities for retaining the measuring accuracies for measurement of structured surfaces. None of the methods show clear advantages over others.
• (3)  
  On measuring the linear patterns, the sequential profiling adaptive sampling always has distinct advantages over the other methods.
• (4)  
  On measuring the tessellations, the three adaptive methods show their advantages on measuring the height-related parameters such as the step height. But they have similar capabilities as other fixed sampling patterns when measuring the lateral parameters such as the mean pitch distance.
• (5)  
  On measuring the rotational symmetric patterns, triangle patch and rectangle patch adaptive subdivision samplings show significant advantages.

Low-discrepancy pattern sampling methods have a similar performance as uniform or jittered uniform sampling. However, it does not mean that they are not an optimizing method; for example, clear advantages have been shown when measuring flatness of flat surfaces [8, 23]. On measuring structured surfaces, their advantages are not apparent and sometimes they may not be a better substitution of uniform methods. Sometimes, uniform sampling may be a better solution compared to other fixed patterns; evidence for this can be seen in table 2, (a).

The fundamental advantages of adaptive sampling are in evidence in this work. These methodologies allocate their sampling efforts according to their earlier sample results or models. In other words, they can adapt the sampling effort to key positions which have higher impact factors on enhancing the reconstruction accuracy than others. Although adaptive sampling approaches have no clear advantages for measuring the pitch distance of crossed-gratings (see table 2, (e)) but have been shown to be effective for other structured surfaces and other parameters of tessellated surfaces.

The challenges, however, of applying adaptive sampling to practical measurements still widely exist. The sequential profiling adaptive method may suffer from the mechanical constraints of stability (e.g. the thermal drift) and accuracy in y-direction scanning. Most of the other efficient sampling methods are difficult to implement within the operation envelope of stylus instruments, with regard to complex scan route designs and redundant scan duration. In terms of interferometers, many of the reviewed sampling methods may be of promise, with the aid of a high-resolution CCD and pixel stratification process or lens auto-switch systems. Considering the positioning errors and optical resolution constraints, a specialized research work on intelligent sampling of the interferometers may be required in the next step. In addition, more theoretical work is necessary to further the research on intelligent sampling. For example, the data storage solutions need to be reconsidered which was introduced at the beginning of this paper. The reverse problem on sampling and reconstruction need to be fully investigated on the basis of geometric measurement. Also, determination of the sample size for an adaptive measurement is a tough research topic which requires particular attention.