Non-equidistant scanning path generation for the evaluation of surface curvature in metrological scanning probe microscopes

In this paper, we present a novel geometry information-based adaptive step (non-equidistance) scanning path generation method for metrological scanning probe microscopes. This method reduces the total amount of required data and enables faster surface scanning speed for large industrial workpieces while preserving adequate geometric information for performance evaluation after surface reconstruction. The grid points are generated iteratively while gaining knowledge of the surface geometry step by step. We focus on the curvature properties and then propose a metric for the curvature information based on the triangulated surface geometry. With certain convergence criteria on the curvature measure variation, the proposed methods promise better surface reconstruction completeness and performance evaluation correctness. Simulations on the algorithm are performed on a typical parametric surface. A brief comparison to height-based scanning algorithm is performed to show the adaptability of the novel method on curvature evaluation. Experimental verifications are conducted to show the efficiency of the proposed algorithm.


Introduction
Quantitative scanning probe microscopes (SPMs), such as the atomic force microscope (AFM), scanning tunneling microscope (STM), and electrostatic force microscope, have been * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. the well-known techniques for surface topography measurements. With the capability of integration with positioning sensors (e.g. interferometers [1]), SPMs are not only widely used as quantitative tools in various fields of scientific research, but they are also developing rapidly for manufacturing and industrial applications, especially for the surface evaluation of microstructures [2][3][4][5][6].
Unlike the attempt to obtain the highest possible resolution in scientific research, the characterization of industrial workpiece samples is conducted over a relatively large region. Further, the efficiency of the measurement either to investigate surface defects or to reconstruct surface geometry has equal importance as the measurement accuracy. However, the measurement speed of large-scale industrial products is often quite slow. On the one hand, such measurements require large datasets for the reconstruction of the measured surface. On the other hand, physical limitations, such as stabilization time for the probe to dwell and read out measurement results, limit the lower bound of the time required for a single readout. A combination of two effects and necessary long total distance movements for measurement ergodicity leads to excessively long measurement times. These long measurement times lead to disturbances on the probe-sample interaction mechanism, making fluctuations in the electronics and variations of the temperature nonnegligible, resulting in unpredictable errors in the measurement results.
One approach to achieve a higher scanning speed is by using non-raster scanning techniques. Such techniques combine circular and linear motions, including spiral [7], Lissajous [8,9], cycloid trajectories [10], concentric circle [11], and optimized Archimedean spiral [12], to avoid driving the actuators with triangular waveforms, providing higher systematic stability to ensure further reduction of the dwelling time. Because of the timespan limitations in the readout electronics and closed-loop controllers, a more commercial approach to reduce the total measurement effort is to sample only necessary grid points. This approach can be achieved by the implementation of adaptive scanning techniques, rather than a conventional raster scan. The selection of measurement grids has been addressed in a variety of studies. Some approaches focus on the precision of the edges [13] using adaptive sampling in successive refinements based on bisection methods. More general approaches reduce the total measured sample points by randomly omitting grid points known as compressive sensing [14,15]. These methods are also often equipped with matching path planning algorithms to offer more suitable scanner movements; thus, the relative positioning errors in terms of the ideal measurement path are reduced [16].
However, current methods mostly focus on height information [13,17,18]. Although it is important for the surface reconstruction algorithm, height information is insufficient for other geometric property evaluation processes in adaptive grid algorithms. Information other than height needs different measures for the convergence criterion. Without such a measure, algorithms tend to oversample the measured surface to extreme, but unnecessary, accuracy in high-frequency regions while ignoring the geometric claims in lower frequency regions. Hence, the convergence criterion considering height information provides a significant refinement effort on steep areas of the surface. For example, in the measurement of Fresnel mirrors [19,20], we show that contemporary algorithms tend to ignore the desired curvature information, which is often in the low varying reflective region, while generating denser grids on sharp edges. Such excessively dense grids challenge both the positioning accuracy of the motion stages and the path planning strategy.
In this paper, we discuss a possible realization of nonequidistant measurement methods for a long-range SPM based on a more practical convergence criterion. To maintain stability during the refinements, we further utilize the curvature measures [21,22] to define the convergence criterion in accordance with the considered geometrical properties (e.g. principal curvatures and normal directions). These convergence criteria are then demonstrated to be more suitable for the maintenance of the local features of the sample surface. Simulations are carried out for a clearer illustration of the proposed method; experimental verification of the method is also addressed by evaluating an optical surface. We show that, with such sampling methods, the reconstructed surface contains sufficient information for the evaluation of its geometric property.

Modeling and methods
The basic concept behind conserving the geometric information of a measured surface is to minimize the difference between the original surface and its approximations, such as triangulation or digitization, while considering certain geometric properties [21][22][23]. These properties, including the area density, mean and Gaussian curvatures, have varying importance in an actual evaluation process depending on various applications of the workpiece. For example, the specific surface area distribution and fractal structure are mostly relevant for biochemical applications [24,25], while surface profiles, such as curvatures and normal vectors, have greater importance in optical microarrays [5,26]. The principal optimization objectives are similar: an adaptive scanning algorithm must include adequate faces to satisfy the performance evaluation of the reconstructed surface for certain applications. In terms of the measurement task, it is equivalent to declare that a sufficient density of grid points must exist to approximate certain curvature properties, and the measured curvature must be stable with respect to the grid refinement.
Conventionally, by choosing a 3 × 3 neighborhood of a grid point and assigning a virtual coordinate system for the grid point coordinates, the surface can be approximated by a second-order interpolation of the measured values on the nine grid points. Let ⃗ z =⃗ z (x, y) be a parametric equation of the measured surface S, and this approximation can be given by z (x, y) = a 00 + a 10 x + a 01 y + a 20 x 2 + a 02 y 2 + a 11 xy, (1) where coefficients a 00 , · · · , a 11 are estimated using the method of least squares based on the measurement readouts of each neighboring grid point [27]. The fundamental forms and curvatures can be derived from this localized parametric surface as an approximation to their real values. However, this approximation does not maintain smoothness between neighboring grids because of the varying estimated parametric surface in equation (1). In particular, this approximation is not suitable when the grids are not equidistant.
Generally, no pointwise convergence of mere principal curvatures exists [28]. Hence, the curvature estimation accuracy of a single point is not strictly dependent on the grid density, especially in terms of the neighborhood of sharp edges. However, the integrals of curvatures of a certain area would, in general, converge to the real value of the sample surface (as (a) Schematic of the curvature measure method. Unlike contemporary methods, the grid point refinement near sharp edges maintains convergence properties, and the tube volume V∆ (B) tends to be its real value asymptotically. (b) Schematic of the tube formula for a specific subset B above surface S. The formula works for ∆ < min ρ, where ρ is the distance between S and the medial axis Sk R d \V to promote one-to-one mapping.
can be observed, for example, using the tube formula [21]), enabling the application of curvature measures as an alternative convergence criterion. Specifically, such curvature measures satisfy at least the O (h ϱ ) , ϱ ∈ R + dependence on gridstep h to promote stability in refined meshes [22]. Such monotone dependency promises theoretically a certain option of grid-step h for arbitrarily high approximation accuracy. The basic concept behind the selection of curvature measures as the convergence criterion is illustrated in figure 1(a), where the normal direction ⃗ u (p) diverges with respect to denser grids, while the nominal volume V ∆ (B) converges to the actual volume above the surface.
However, for unknown a priori surface structures, determining exactly how well an approximation conserves the necessary information is often difficult. One possible approach is to develop a similar convergence criterion to the pure height measurements presented in [13,18]. This is based on using the measurements from two prior refinement levels to verify whether the surface converged locally. Wherever the difference between the estimated surface curvature measured between two prior refinement levels is larger than the threshold, the scan path is refined further based on the convergence dependence of O (h ϱ ). In this section, we state the main ideas of the approach and focus on the basics and derivations. Then, the refinement based on the O (h ϱ ) dependency and the Laplacians on triangle meshes [29] is modified accordingly.

Notation and key concepts
Before addressing the convergence criterion in the measurement process, a series of fundamental concepts that will emerge frequently in the following discussion should be presented.
A surface is denoted by S in the three-dimensional (3D)oriented Euclidean space R 3 . A surface patch σ (u, v) is a localized parametrization of the surface with (u, v) ∈ R 2 . The values of the unit normal vectors ⃗ n at the points of S (or σ) are recorded by its Gauss map G s , which is the map from S ∈ R 3 to the unit sphere S 2 that assigns to any point p ∈ S a unit vector ⃗ n p ∈ S 2 . The rate at which ⃗ n varies across S is measured by the derivative of G which is often known as the Weingarten map. The Weingarten map can be symmetric; the associated quadratic form is referred to as the second fundamental form. The eigenvectors and eigenvalues of the Weingarten map are called principal directions and principal curvatures, respectively. Both principal curvatures, denoted by κ 1 and κ 2 , can be recovered from the trace and determinant of D p G. The mean and Gaussian curvature at p are given by 1 2 (κ 1 + κ 2 ) and κ 1 κ 2 , respectively.
For closed-loop surface metrology systems, such as constant signal mode STMs and AFMs or micro-coordinate measurement machines, the probe center is constrained to a surface of constant separation away from the original workpiece. Let ∆ be some unchanged separation corresponding to this scheme, and the parallel surface S ∆ of S is defined as [30] S ∆ = {p + ∆⃗ n p |p ∈ S}. (3) Let V be the workpiece in R d , and then S denotes the surface of V to be measured. Let ρ be the distance between S and the medial axis Sk of the complement of V in R d and V ∆ denote the volume between S and S ∆ . Assume that ∆ < ρ for all points p ∈ S (the parallel surface S ∆ does not intersect with itself such that the probe is assumed to measure only one point once). Then, for a subset B on S, we obtain the tube formula [31] Vol where H and G are the mean and Gaussian curvatures, respectively, and H 2 is the 2D Hausdorff measure. Figure 1(b) illustrates a schematic of the tube formula. This inspires the introduction of curvature measures as a measure of the polynomial dependence on ∆ of the tube volume of a specific subset B ⊂ ∂V = S. The concept of curvature measure [21] is a generalization of the Gaussian and mean curvatures for convex and smooth objects. It was subsequently extended to triangulations, digitized objects, and sub-analytic sets by P. Wintgen [32] and M. Zähle [33] by the introduction of normal cycles, which is the integration current defined by the graph of the Gaussian map in equation (2), and it serves as basics of the generalization of the curvatures in non-smooth geometry [34,35]. The curvature measures are typically denoted by µ k , and in the smooth case, we have µ 0 = dH 2 , µ 1 = −2HdH 2 , µ 2 = GdH 2 and µ 3 = 0.
Specifically, to ensure a clear description of the following discretization method, we follow the subscript conventions given by [22,36], rather than by [34], where the subscripts are in reversed order. Therefore, the volume given by equation (4) for B in 3D manifold can be rewritten as where µ k (B) denotes the integral of µ k over B and ω p denotes the volume of the unit p-dimensional Euclidean ball. Below, we list the values of ω p for small p.
Thus, we are able to assign µ k (B) to each individual subset of the surface B ∩ S. To recover the classical geometric properties, the curvature estimators are also defined. The normalized mean (Gaussian) curvature estimatorĤ u (Ĝ u ) of S is given bŷ The normal vector field u can be represented in various forms based on different vertex normal computing methods. A possible method is the Voronoi Covariance measure [37,38], which provides a convergent normal vector field ⃗ u for a certain family of parameters. A carefully chosen approximation method of the normal bundle can improve the estimate significantly, leading to more accurate curvature measure evaluation. In this study, the general polygonal meshes during the refinement were triangulated implicitly by splitting each nontriangular face into triangles at its barycenter for simplicity, and the normal vectors on the newly added grid points are estimated by interpolation of the nearby facets. As highlighted in [22], such simplifications are not influential to the converging stability. The class of such normal vector estimator selections are called the corrected normals. Thus, the corrected normal vector evaluations yield the corrected normal cycle and, consequently, the corrected curvature measures.
After specifying the vertex normal, the closed form formulas for the curvature measures can be derived easily. (x i ) i=1...n denotes the positions of the n grid points of the surface polygonal mesh. Let (u i ) i=1...n be some normal vector field prescribed at these points. Therefore, on an arbitrary triangle τ ijk on the generated mesh with vertices i, j, and k, the interpolated corrected curvature measures take the following forms: where ⟨·|·⟩ denotes the scalar product andū = 1 3 (⃗ u i +⃗ u j +⃗ u k ). The measure µ (0) is the corrected area density of the given triangle, and µ (1) and µ (2) are the corrected mean and Gaussian curvature density, respectively. Further, µ (X,Y) is the trace of the corrected second fundamental form along arbitrary directions X and Y, as an alternation of µ 3 in equation (5) for corrected normal ⃗ u. µ (X,Y) can be used to determine the distance of the corrected normal ⃗ u from the ground normal ⃗ n, as a measure of the quality of the choice of corrected normal. When ⃗ u = ⃗ n and h → 0, µ (k) (B) → µ k (B) , k = 0, 1, 2 and µ (X,Y) (B) → µ 3 (B) = 0 is obtained asymptotically. The ground normal ⃗ n requires complete knowledge of the surface profile and is unrealistic for measurement tasks. However, it can be used in a simulation processes as a primary verification of the algorithm.

Adaptive measurement algorithm
The main concept behind this approach is to minimize l 2 errors of the estimated curvature measures between two consecutive measurements. The key points of the refinement algorithm are as follows: (a) Perform measurement on a net formed by rows and columns with predefined h discretization. Then, estimate the concerned curvature measures pointwise using the aforementioned curvature measure method. The estimated results are denoted by Meas  (k) ii of the same dimension, whose differences can be calculated pointwise. Identify grid points where the estimated curvature measures between the last two iterations differ more than a given l 2 error criterion. (d) At regions where the convergent criterion fails, the net is refined twice in resolution, to h 2 . The algorithm repeats from the first step until the desired convergent criterion is satisfied everywhere. Figure 2(a) depicts a schematic illustration of this h 2 estimation process. Here, the black dots are grid points of the first performed measurement, which are then interpolated to h 2 separation, as depicted by the dashed line at the center. Another measurement with h 2 displacement in both X and Y direction was performed in the next iteration and represented by blue dots. The actual value in the second measurement (blue vector) is compared to the estimated value in the first measurement (gray vector), and the error is the difference between the vectors. In a triangulated setup, the h/2 grid deviation process designed for rectangles is more practical to be replaced by deviating to the barycenters of the triangle faces to avoid selfedges and repeat face elements. Notably, the measurement path must be generated to follow each refinement step. Potential candidates for such paths are the space-filling curves with a fractal dimension of two [39], considering the bisection nature in subdivision spaces. A representative of such curves are the 2D Hilbert curve and Sierpiński curves, whose quadtrees subdivide the side length of the squares into two parts in each step [40]. The triangulation of such distributed sample points is also well-studied [41], and the evaluation of the curvature measures can be given directly from the triangulation result. Moreover, a routing system based on the space filling curves is a promising solution for the planar traveling salesman problem (TSP) of distributed individual points in the scanning path by equalizing the weights of all refined points [42,43].
We can derive the interested surface properties from the previous defined curvature measures by choosing the integration base B ∈ R 3 . One possible approach is to apply the initial grid as the integral limits (i.e. a square with a side length of h and extruded by ∆ on both sides, as shown in figure 2(a)). A more flexible approach introduced in [22] is to define a ball around⃗ x, which may or may not be larger than the initial grids, but covers the sharp edges better and significantly cleaner in mathematical expressions. Let B r (x) be the ball of center ⃗ x and radius r. The corrected mean curvatureĤ and the corrected Gaussian curvatureĜ at ⃗ x are defined in accordance with equation (7) byĤ The first and second principal curvaturesκ 1 andκ 2 at ⃗ x and their associated corrected principal directionsv 1 andv 2 are also defined bŷ where are the eigenvalues of M, and ⃗ z 1 , ⃗ z 2 , and ⃗ z 3 are their associated eigenvectors.
is the Kronecker product of the corrected normal ⃗ u with itself, and some large coefficient K is chosen to force the tangency of principal direction eigenvectors. We are also able to implement alongside another iterative refinement algorithm considering the determination of sharp edges. The corrected normal cone algorithm is also well suited for the convergent detection of such discontinuities, as shown in the schematic illustration in figure 1(a). The key points of this complementary algorithm can be given as follows: (a) Define sharp edges and extreme points as regions where the estimated curvature exhibits distinct statistical properties, compared with the majority of the sample surface, whereas the curvature estimation follows the same steps as given in the above refinement algorithm. (b) Novel refinement points are added to the barycenter of triangulations with distinguished properties, and the integration ball radius follows the same bisection asymptotical reduction criterion as that given in the above refinement algorithm. Add these points in the next curvature value estimation process, and if these points are indeed extreme points, the total area of such distinct distributions is supposed to be suppressed face-wise. (c) When the total area of the edges and extreme points does not shrink further, the algorithm reaches a convergence criterion, either on the difference between two consecutive measurements or a total portion with respect to initial triangulation.
Due to the implementation of this second refinement algorithm, the system can distinguish the edges from the surfaces. Further, the surface to be estimated is no longer limited to a predefined region, but expands along with the shrinkage of the edge area in the successive iterations, leading to more accurate curvature estimations near the neighborhood of these extreme points. The integration ball r can also be refined at each iteration step to discover smaller features in accordance with the finer grids. A flowchart for the combination of the aforementioned algorithms is shown in figure 2(b).
As a conclusion, the total measurement time can be formulated approximately as a combination of the measurement time of the sample points Nt a , the traveling time along the generated paths (N − 1) t s and the calculation time for obtaining the geometric information T c , given by: where the footnote i represents the ith iteration, N denotes the number of the generated refinement points, t a and t s represents the average data acquisition time and the average traveling time between consecutive points, respectively. The proposed algorithm contributes to saving the total measurement time from all the three aspects:  [44]; so the total calculation time for retrieving geometric information is also reduced due to a smaller N.
We exploit this combined non-equidistant grid point generation algorithm to perform adaptive measurement in the following sections.

Numerical simulation
In this section, we implement the aforementioned combined refinement algorithm into the virtual measurement process of a multi-grade parabolic Fresnel mirror surface profile to illustrate the adaptability and convergent performance of our geometric property-oriented adaptive refinement algorithm. We first provide a primitive setup to show the adaptivity of our method on the evaluation of curvature properties. Then, a more detailed simulation process is presented to more clearly illustrate how the presented algorithm works in iterative refinement steps.

Numerical simulation methods.
Following the steps given in the previous section, a first step refinement of the grid points of our algorithm and the contemporary height convergent method are both shown in figure 3 to provide an initial perception of our method in contradiction to the contemporary height convergence method. The initial data points are generated in figure 3(a), with a separation of 0.01 unit length. The two measurement nets of the first iteration are then generated with mesh grid of 0.2 unit length, with a displacement of 0.1 unit in both directions, as depicted in figure 3(d). The     (1) ii , the conclusion is straightforward: the height error convergent criterion, which is not explicitly responsible for curvature profiles, are not sensitive to the curvature errors and will not consider these errors during further refinement steps. Meanwhile the curvature error convergent criterion addresses such errors properly, as denoted by the red triangles in figure 3(f). This shows the adaptability of the modified algorithm on the evaluation of curvature properties.
A more detailed simulation process is presented for the combined algorithm, where the iterative refinement and the number of the data points are also considered. The surface is triangulated for simplicity in the integration process, and the l 2 convergent criterion is replaced by l ∞ error, which is the maximum deviation of the estimated curvatures in consecutive measurements. Determining the edges leads to a shrinking portion of the area of the edges and an expanding portion of the area of the surface patches in the total projected area on the uv-plane. Figure 4 shows the simulated first-step refinement setup of a three-level Fresnel mirror. Grid point separations of 0.4 and 0.2 mm are generated. Figures 4(a) and (b) provide the corresponding estimation of the mean curvatures. By fitting the estimated curvature to the theoretical surface, which is often known for manufactured species, or roughly assume that the majority of the estimated values must follow a normal distribution profile, we can regard those extreme values as initial guess for the position of the edges. Consequently, by setting the edge determination threshold to 3σ with σ the standard deviation of the fitted normal distribution, figures 4(c) and (d) show the surface meshes after removing the identified edges from the initial triangulation. A comparison of the curvature is performed in accordance with figure 2(a), where the estimated values Meas (1) ii of the refined points at the barycenters (denoted by the red dots in the subfigure) of the initial triangulations are interpolated by an integration within the integration ball of radius r, as also shown in figure 4(c). After the refined points are measured, a novel triangulation, including these barycenters and the refined edge areas, is generated, and the corresponding Meas (1) ii is evaluated with the same integration ball. After a comparison of the two datasets Meas (1) ii and Meas (1) ii , barycenters with deviations larger than a certain threshold are retained in the next iteration. The edge regions are herewith re-evaluated, and the regions marked as the surface profile are also re-generated for the evaluation in the next iteration. The integration radius is decreased iteratively to capture smaller and smaller surface features. The refinement processes of both grid separations are shown in figure 5, whose subfigures (a) and (b) show the first refined grid and the refined grid after four iterations for grid separation of 0.4 mm. Figure 5(c) illustrates the estimated mean curvature distribution of the final refinement, in comparison with the estimation results given by figure 4(a). Similarly, the (d) and (e) show the grids for 0.2 mm, and (f) depicts the corresponding estimation of the mean curvature across the entire surface. The convergent performances of both grids are shown in figure 6.
Another example is simulated on the surface of a multifocal micro lens array (MFMLA) [45]. The results are given in figure 7. The subfigures (a) and (b) show the refinement process of the proposed curvature convergent algorithm and illustrate the surface profile of the MFMLA used in the simulation. Same as the previous examples, the measurement process

Simulation results and discussion.
The simulation results indicate the convergence property of the iterative grid point generation algorithm. Figure 4(a) shows the 0.4 mm grid setup hardly conveys accurate information about the mean curvature properties. This is also partially true for the 0.2 mm grid, where, although the initial estimation of mean curvature  can distinguish the boundary of the fast-variating regions from smooth surface patches, it fails to resolve smaller surface patch profiles in between, as expected. This leads to overly large vacancies of the identified edge regions. However, after four iterations of grid refinement, the edge information is already able to be derived from the estimated mean curvature, as shown in figure 5(c). When the initial grid separation is set to a smaller value (0.2 mm), the final estimation of the mean curvature addresses the edges better. Zoomed images of the final refined grids of the two different initial grid separations in figures 5(b) and (e) show that the coarser grid spares more effort on the refinement of grids near sharp edges in successive iterations, while finer initial grids indicate a much-relieved edge refinement strategy. Figure 6 shows the convergence performance of the two refinement processes. Notably, in the first iteration of the coarser grid setup, the estimation of the mean curvature, as well as the determination of the edges, are less accurate, resulting in unusually suppressed data at the beginning of the line chart. The later iterations are relatively free of this type of unconformity, and the convergence nature of the combined algorithm is illustrated as the edges are addressed more accurately in refined grids and the integration radius is set lower to capture more detailed variations. Although an absolute convergent criterion, such as a threshold on the mean curvature deviation or the edge area portion, can be addressed directly in accordance with the simulated convergence plot, it is still strongly recommended that a relative convergence criterion focusing on the difference between consecutive measurements be used instead, especially with unknown manufacturing quality and systematic measurement errors.
Further, the edge determination threshold strongly affects the refinement process. Table 1 shows the influence of different criteria used when determining the extreme points that are supposed to be edges. For a looser threshold of 3σ, which is used for the illustrative simulation results, the number of  Threshold  Iteration  1  2  3  4   1σ  Surface point # a  289  501  786  1173  Edge point #  788  3328  10597  27664  3σ  Surface point #  892  757  1250  3465  Edge point #  324  2919  8504  15700 a The numbers represents only the newly added number of grids in the corresponding refinement iteration; the total number of measurement points is a sum of the rows and the generated initial grids.
the newly added face refinement points grows more rapidly in comparison with a stricter threshold of 1σ. However, the number of the refined points in the edges with 1σ threshold grows even more drastically. This is partially due to a larger number of triangulation elements that are marked as edge elements. Another reason is due to the surface model itself: by theoretical calculation, a significant portion of surface curvature lies near or beyond the [−σ, σ] interval centered at the universal mean. This indicates the importance of the knowledge on the measured surface model for better time-saving convergence criterion choice strategy. Now consider the simulated refinement process of the MFMLA as shown in figure 7. The applied height convergent algorithm also fails to capture detailed curvature information, as shown in the right-hand-side of figure 7(d), especially within regions neighboring to flat surface patches; the proposed curvature convergent algorithm, on the other hand, is able to capture such discontinuities properly. But there are also issues on the curvature calculation process: the computation on rims and corners of the measurement region may fail due to the application of smaller integration ball radius, as shown in the zoom-in illustration. This is not a problem when the measurement scanning range is chosen to be larger than the feature region; but for totally unknown surface profiles, the integration ball radius should be chosen larger to overcome such failures. The total number of measurement points in the proposed algorithm is 48 134, which is also larger than a total of 39 374 points in the contemporary height convergent algorithm. This problem can be partly addressed by selecting weaker criterions, but the accuracy on the curvature estimation also becomes worse.
A comparison of the convergent characteristics of both algorithms with different measurement concerns is also given in figure 7(e). As can be seen from the results, although the proposed algorithm outperforms the height convergent algorithm in the curvature estimation applications, its performance on pure height measurements is less effective, as expected.
In general, the above simulation process shows the theoretical convergence performance of the proposed algorithm. The algorithm shows adaptability in the evaluation and estimation of surface curvature, compared to the height convergent method, which often ignores the diversity in such surface properties. The algorithm also shows relative robustness in that the dataset used contains artificial noise. This robustness benefits from the integration nature of the corrected normal cone method. Codes used in the above simulation process can be found in the supplementary materials (available online at stacks.iop.org/MST/32/125009/mmedia) as well as our Git-Hub repository: https://github.com/yyhuzju/non-equidistantalgorithm/tree/no_GUI.

Experiment setup
In a previous study, we developed a metrological STM for the measurement of large area workpiece measurement [3,6]. A primitive verification of the proposed algorithm was performed on this STM system. A schematic illustration of the system setup is depicted in figure 8(a), and a photograph of the system is shown in figure 8(b), illustrating the relative assembly relationship between the actuation platforms. The STM probe was prepared via electrochemical etching process [46]. Figure 7(c) shows a photograph of the probe tip acquired using scanning electron microscope. The probe has a sharp tip with a radius below 50 nm. The adapted linear positioning stages (Type M-112.1DG1, PI Inc.) for the x-, y-, and z-axis have a minimum incremental motion of 0.05 µm and a unidirectional repeatability of 0.25 µm, with a traveling range of 25 mm. The position of the scanner axes during the measurement process is acquired by built-in encoders. A two-arm 632.8 nm wavelength laser interferometer was used to measure the translational volumetric errors and the errors are compensated according to a geometric error model. The STM probe was installed on a piezoelectric actuator (PEA, Type VS15, CoreMorrow) with a traveling range of 19 µm. All the stages were mounted on a marble bed. The tunneling current was measured using a picoampere electrometer (Type 6514, Keithley Instruments), and the readouts were fed directly to a PC runtime controller via the IEEE 488.2 protocol.
The dwelling time of the probe is given an upper limit for a single readout to converge, which is 2 ms for points in the identified regular areas and 5 ms for points on edges. Readouts that are not convergent within the given time limit are approximated by averaging the last ten results. This kind of rough cutoffs are usually less influential in the case of curvature measurements than in height measurements, because of the integral nature of the curvature estimation method.
The workpiece to be measured is a piece of Fresnel reflector array with a parabolic profile. Figure 9(a) demonstrates how the Fresnel reflector array works: for parallel incident light, the reflector focuses the light on a plane above the reflector surface to distributed focal points. The images were acquired   by an optical microscope and showed both the sample surface and the focus plane, whereas the later one was actually distributed images of the light source. The cross-section profile is also given in figure 9(a), showing the dimensions of the multistage parabola. An on-machine setup of the tip-sample pair is shown in figure 9(b). When performing the measurement, the separation between the probe tip and the sample was set to approximate 0.05 µm to avoid collision and damage to the workpiece surface.
The measurement path follows the outstanding spacefilling curve solution of the TSP [42]. The reposition point was set to approximate the origin of the workpiece. A primitive comparison between the total scan distance is shown in figure 10. The figure indicates that, by using the space-filling curve solution, the total traveling range of the scanning probe is reduced to 63.49 mm with 31.21 and 32.28 mm for the original and the deviated grids, respectively. Meanwhile, in the traditional raster scanning scheme, the total tour range is 115.30 mm with 56.58 and 58.72 mm for the same two grids. This shows the efficiency of the space-filling curve solution in reducing the total movement time for faster data acquisition. Further, space-filling curves promote more tender corners, thus avoiding instability in fast positioning during the measurement process. A more general solution is by 3D space-filling curve with knowledge of the surface profile from former measurements. We chose 2D curves in this study for simplicity, as the vertical height difference is relatively small compared to the long horizontal traveling range.

Measurement results and discussion
The measurement results are shown in figure 11. Figure 11(a) shows the result under the finest grid that can be measured using the current apparatus. The grid separation was set to 5 µm, and the outermost points of the cylindrical workpiece are extrapolated to form a square area for better curvature estimation purpose. The mean curvature is estimated in accordance. A zoomed-in scope showed the processing traces by the diamond tool. The initial estimation with a grid separation of 0.2 mm as a start point of the refinement process is shown in figure 11(b). After six iterations of refinements, the estimated mean curvature of the surface is shown in figure 11(c). The 3D model of the final measurement is also given in figure 11(d).
In general, the algorithm used in this study fits well for the estimation for the curvature estimation process. The subfigures show that the estimation result from the first measurement is blurred, and thus it is difficult to obtain useful information. However, after several iterations, the major surface features are clearly visible and a curvature estimate is obtained similar to that of the densest mesh. However, the processing traces are also neglected due to the coarse grids in the flat regions. Finer initial grids may help identify such small features, but may also lead to longer acquisition time. The number of points added in each iteration is given in figure 12(a), which is 138 944 points in total. For comparison, the total number of points in the finest raster grid is 246 253. This reduction in the total measurement points yields an over 30% improvement in the total measurement time with space-filling curve paths, and it is expected to increase the measurement efficiency by over 50% for the raster scanning scheme. The total approximated measurement time is given in table 2 as a reference. Note that in more advanced setups, the total traveling time can be reduced further by applying faster scanners (i.e. reducing the average t s in equation (13)), then the upper limit of the total measurement time will be governed by the data acquisition (setting) time Nt a . In sparser grids, however, such as those in the first few refinement iterations, the movement of the scanner still needs time quite longer than the data acquisition process, i.e. t s ≈ 150 ms for average 0.2 mm point separation when the acquisition time remains the same as given above. In this situation, the path planning is of the same importance as the adaptive grid selection. As a result, the total traveling time after path planning is ≈ 2393 s, while for raster grid with the shortest average traveling time of ≈ 0.012 s for finest grids as given in the table, this  time will be ≈ 2955 s. The calculation time in the table refers to two processes: solving the TSP and computing the corrected curvature. The codes given in this paper is not optimized and the integration process is rather time-consuming. This problem can be partly addressed considering that when retrieving curvature information, the integration process is unavoidable, and the calculation time of 138 944 points implemented in the iterative process is still faster than the same process on 246 253 points after performing traditional raster measurements. Thus, despite some estimation noise in the measurements, the algorithm proposed in this study still has a considerable practical value in the curvature estimation of complex workpiece surfaces. The path connecting the refinement grids in the first few iterations are shown in figure 12(b). The path in each iteration starts from the origin, travels toward the first quadrant, follows a clockwise direction, and finally ends at the same origin, for the convenience to wait for subsequent instructions. The minimal point separation in the 0.2 mm case is approximately 3 µm, which is still acceptable regarding the specifications of scanner. The mean curvature estimation uncertainty for consecutive five measurements is given in figure 13, showing the histogram of the estimated mean curvature by area sum. As can be seen from the results, the estimation in the regular regions, i.e. the total area with smaller mean curvature, is relatively stable for different measurements. However, the side bands around the central peak in the finest grid are suppressed. This is mainly because that the small variations in the flat areas such as the processing traces by the diamond tool shown in figure 11(a) are averaged and neglected in the adaptive grid. This problem can be solved easily with finer initial gird covering the feature size of such traces, but for large area profiling, such problems have limited impact.

Conclusion
In this paper, we proposed a novel convergent criterion on the adaptive grid generation algorithm for SPMs. For surfaces with large regular areas, the use of this adaptive scanning algorithm, rather than using traditional equidistant grid scanning, can save the total measurement time without losing significant information about the properties of interest. Furthermore, unlike the contemporary height convergence algorithm, this method captures the curvature information and controls the estimation accuracy of surface profiles related to curvature properties. Further, the combined algorithm with edge identification can distinguish edges from regular surfaces, leading to more accurate identification of fast surface variations. The contribution of each part of the algorithm can be termed as follows: • The iterative point refinement process saves the total traveling and data acquisition time simultaneously by reducing the total number of points to be measured. • The TSP solution used in the path planning process can reduce the traveling time further, especially when the surface feature points are sparse and the total movement time between consecutive sample points governs.
We also presented a simulation process to theoretically show the convergent performance of the algorithm and verified the algorithm on a micro-structured surface by using an iterative measurement process with exceedingly higher resolution on the critical surface regions, performed on a self-built metrological STM. The results demonstrate the efficiency of the developed algorithm in reduction of the total acquisition time. Equipped with the space-filling curve solution for the shortest scanning path, the total measurement time can be reduced further, and the stability of the entire system is improved.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.