The measurement of texture for geometric surfaces is well established for surfaces that are of a planar (Euclidean) nature. Gaussian filtering is the fundamental base for scale-limited surfaces used in surface texture, but cannot be applied to non-Euclidean surfaces without distortion of the results. A link exists between Gaussian filtering and solutions of the PDE that models linear isotropic diffusion. In particular, an analytical solution of this diffusion equation over a planar region at a time t is given by the continuous convolution of the initial distribution of the diffused quantity with a Gaussian function of standard deviation . A practical implementation of the standard Gaussian filter on sampled data can be viewed as a discretization of this process. On a non-Euclidean surface, the diffusion equation is formulated by using the Laplace–Beltrami operator. Using this generalization, a method of Gaussian filtering for freeform surface data is proposed by solving the diffusion equation for approximation residuals defined on a freeform least-squares approximation of the measurement surface data. Results of the application of these methods to simulated and experimental data are presented.
Surfaces provide the functional interface through which all advanced and emerging products operate at the most fundamental level. In all cases, functional performance depends on our ability to assess the quality of the surface topography/geometry.
During the last decade, measurement of texture for geometrical surfaces is mature for ‘planar’ (Euclidean) surfaces: that is to say for surfaces that can be defined by a single height value for each point in a plane. Characterization of surface texture over a planar area is well established (Jiang et al. 2007a,b; ISO/FDIS 25178-2 2010). Filtering techniques for ‘planar’ (Euclidean) surfaces based on Gaussian, spline and wavelet techniques have been comprehensively investigated during the last decade (Jiang et al. 2007a,b; ISO/TS 16610 Series 2010). Gaussian filtering is the standardized base for scale-limited surfaces used in surface texture (ISO 11562 1996; ISO/CD 16610-61 2010; ISO/FDIS 16610-21 Series 2010; ISO/FDIS 25178-2 2010).
However, with the development of science and technology, the current industry has grown from one based on skills and manual labour to one based on design and manufacture of advanced products. For example, in the optics industry, the high-added-value product spectrum has shifted to design and fabrication of freeform optical elements and micro-structured optical surfaces that are critical components in photonics, telecommunications, laser printers, scanners, displays and broadband optical fibre connectors. Unlike conventional surfaces, the surfaces of these advanced components have no axes of rotation and no translational symmetry and in the future could have almost any designed shape. Surfaces that have no rotational or translational symmetry have been termed freeform surfaces (Clayor et al. 2004; Jiang et al. 2007b); in particular, they mostly require sub-micrometre form accuracy and sub-nanometre surface texture.
Characterization of surface texture on such complex freeform geometries at micro-nanoscale is proving very challenging. Apart from the problem of best fitting the nominal geometry to the measured surface and subsequent calculation of the residual surface, there is a problem of how to decompose the residual surface into different scale-limited surface components (if we take a traditional approach, this means how to separate residual surface into roughness and waviness components).
In most real-world applications, including surface metrology, surfaces take the form of a two-dimensional manifold embedded in the Euclidean space . When the measurement surface in question is a plane, this is still the case; however, it is a specific case of a two-dimensional manifold having zero curvature everywhere and is therefore also Euclidean. For non-planar surfaces, the manifold contains points that have non-zero curvature, and it is therefore a non-Euclidean space. In the Euclidean case, given a line and a point, there exists only one line through that point such that it is parallel to the first line. This is known as the parallel postulate and is one of the many equivalent definitions of Euclid’s fifth postulate. However, for a freeform surface measured by a lattice method, there can be many (or no) geodesics through a point that are parallel to another (i.e. non-Euclidean geometry). This property is illustrated in figure 1. However, a freeform surface cannot transform to a plane that retains both orthogonality and scale. Therefore, current surface characterization techniques (defined on Euclidean geometries) are distorted. This can be illustrated with the following example. Consider a hemispherical surface where measurements have been taken on a grid relative to the instrument plane. The surface texture can be represented by the deviations between the measurement data itself and a hemispherical surface that best approximates the data. At points near the hemisphere boundary, we have a case where the tangent plane of the hemisphere is almost perpendicular to the instrument plane. If the surface texture were to be displayed on the measurement grid, the flattening process will cause compression of the surface and will result in smaller implied distances between neighbouring points than their actual geodesic distance on the surface. Thus, this planar projection of the texture is distorted as would be any Gaussian filter cut-off value and surface roughness parameter obtained from it. Also, this distortion would be different at different points on the surface. Furthermore, the planar projection method breaks down when the surface starts to go back on itself as in all closed surfaces.
In this paper, we provide an initial investigation, which uses a specific PDE (that has solutions theoretically linked to linear Gaussian filtering) to filter data. In Euclidean geometries, this method is equivalent to linear Gaussian filtering, but it is possible to use the theory of differential geometry to solve the PDE on non-Euclidean surfaces. This theoretically provides a tool for performing ‘Gaussian filtering’ on a surface that is appropriate for and dependent on the geometry of that surface and should yield results that are free of distortion. The nomenclature used in this paper is given in table 1.
2. Euclidean surfaces
In surface metrology, nominal surfaces are regarded mathematically as a continuous function that defines a height value over a planar domain of interest. Formally, if the domain of interest is denoted by , then the surface height is given by the function , and the continuous surface is described by the graph of g which is the set2.1As the domain of g is Euclidean, well-established methods for analysing the surface data are available (Fourier analysis, wavelet decomposition, linear filtering). The Euclidean nature of such surfaces means that these analytical techniques can be applied to the data without need for theoretical consideration of the effects of curvature.
3. Freeform surfaces
Many types of simple surfaces (such as planes, spheres or cylinders) can be described by a simple equation. For more complex surfaces, such as those without rotational and translational symmetry, it may not be possible to obtain a simple formula to describe them adequately. Such types of surface are referred to as freeform. Although it is possible to obtain an approximate functional representation of freeform surfaces by fitting approximations to a set of discrete measurements taken from the surface, from a point of freeform design and manufacture view, the popular choices of approximation form for this purpose are parametric surfaces such as B-spline and non-uniform rational B-Spline (NURBS) surfaces (Piegl & Tiller 1997; Jiang et al. 2010).
As parametric surfaces are referred to a great deal in later sections, some notation is required. A general parametric surface is denoted by and is parametrized by a function such that3.1where is the parameter space. Up to now, the filtration of surface metrology data has been performed using a variety of techniques including wavelets and Gaussian filters. These techniques assume that the underlying geometry of the surface portion being analysed is planar and therefore exhibits no curvature. In the case of measurements from freeform surfaces, this assumption is no longer valid, and applying the same filters to such data will result in distortion as a result of the curvature. In particular, band-pass filters based on Fourier analysis that are valid for Euclidean geometries are no longer immediately applicable on geometries that are not Euclidean. If progress is to be made in the area of filtering freeform surface measurement data, techniques that yield reliable and valid results for non-Euclidean geometries are required.
4. Linear Gaussian filters
The linear areal Gaussian filter is the standard method for analysing measurement data of surface texture for nominal planar surfaces. Linear areal filters have an associated weighting function s(x,y), which in the case of the areal Gaussian filter is given by4.1where x and y are the distances from the centre (maximum) of the weighting function in the x- and y-directions respectively. λc is the cut-off wavelength, and α is given by4.2The value of alpha is chosen to provide 50 per cent transmission characteristic at the cut-off wavelength λc. The filter is applied by performing a discrete convolution with the surface data and discrete values sampled from the weight function s(x,y).
Repeated application of the Gaussian filter using an increasing sequence of standard deviation parameters yields a set of increasingly smooth representations of the initial data, known as the Gaussian scale space (Witkin 1983). This provides a means of analysing data at a number of increasingly coarse resolutions.
5. Diffusion filtering and relationship with Gaussian filters
(a) The diffusion equation
Diffusion equations are used to model changes in concentration of a quantity of interest inside a specified region with respect to spatial and temporal variables. If the value of the quantity of interest at a particular point in space and time is described by a continuous function (where denotes the volume of interest), then the evolution of this quantity is described by the PDE5.1where Δ denotes the Laplacian operator with respect to the spatial variables. Equation (5.1) is more commonly known in the PDE literature as the heat equation and is a special case of the more general anisotropic diffusion equation (Weickert 1998). As mentioned in §2, a general continuous surface defined over a planar domain can be represented by equation (2.1). This is the graph of a function where .
If a boundary condition of the form f(x,y,0)=g(x,y) is imposed on the diffusion equation (5.1), it is known that its solution f(x,y,t) at time t is given by a continuous convolution of the function g(x,y) with a Gaussian function of standard deviation . The analytical solution of the diffusion equation is expressed as a continuous convolution, but in general practice, a measured surface is in a discrete format and the function values to be diffused are discretely sampled. In the discrete setting, a common approach is to approximate the convolution process by performing a discrete convolution using sampled values from the Gaussian kernel (Weickert 1998). Furthermore, performing diffusion filtering on a freeform surface can be achieved by generalizing equation (5.1) to non-Euclidean manifolds with the use of the Laplace–Beltrami operator (Desbrun et al. 1999; Clarenz et al. 2000; Bajaj & Xu 2003).
(b) The relationship between the diffusion time parameter and the Gaussian cut-off wavelength
The relationship between the time parameter t of the diffusion process and the standard deviation σ of the equivalent Gaussian function is given by . From this relationship and the definition of the areal Gaussian weighting function (4.1), the required value of σ can be calculated for a given cut-off wavelength λc according to5.2From this relationship, the appropriate diffusion time parameter t is given by5.3In the Euclidean case, the cut-off wavelength provides a 50 per cent transmission characteristic of a sine wave. In the case of freeform surfaces, since there is currently no clear definition of the equivalent to a sine wave on a freeform surface, this interpretation requires further exploratory research. It is sensible though to interpret the value of λc as a nesting index for data smoothing until a fully developed interpretation is available.
6. The Laplace–Beltrami operator
The Laplace–Beltrami operator is the generalization of the Laplacian operator to functions defined on surfaces or more generally Riemannian manifolds. When the manifold in question is Euclidean space, the Laplace–Beltrami operator simplifies to the standard Laplacian operator. Although the definitions apply to general manifolds, for areal surface filtering, it is two-dimensional manifolds embedded in that are of interest, particularly parametric surfaces such as NURBS surfaces. In order to define the Laplace-Beltrami operator, some basic definitions from differential geometry are required. To allow a more compact description, partial derivatives are denoted by6.1where i=1,2; j=1,2; u1=u and u2=v. Then, using the notation introduced in equation (3.1), consider a parametric surface and a function such that . The tangent vectors at a point are given by xu1 and xu2, from which the surface normal n at p can be evaluated as6.2using the Euclidean norm ∥.∥ and vector product ×. The first fundamental form of is denoted by the matrix and its components are given by6.3where i,j=1,2, and the inner product 〈·,·〉 is the usual Euclidean scalar product. We can now define the tangential gradient operator acting on a function by6.4and the tangential divergence operator acting on a C1 continuous vector field r on by6.5The Laplacian operator on denoted by is then defined as (do Carmo 1993)6.6which can be represented explicitly in terms of coordinates by6.7With equation (6.6), we can now model the isotropic diffusion of a quantity described by a function f defined over the surface as6.8
7. Mean curvature motion
For any point on the surface, it is well known from the theory of differential geometry (do Carmo 1993) that the Laplace–Beltrami operator is related to the mean curvature according to7.1where H(p) and n are, respectively, the mean curvature and surface normal vector at p. This can also be thought of as applying the Laplace–Beltrami operator to the identity function on the surface . Substituting this relationship into the diffusion equation (5.1) then leads to7.2which is commonly referred to as mean curvature motion or geometric diffusion. Equation (7.2) essentially describes evolution of a continuous surface with a velocity normal to the surface at each point and having magnitude equal to the mean curvature. In the discrete setting where data take the form of a mesh, solving the diffusion equation at increasing time steps will result in a sequence of increasingly smoothed surface meshes. Thus, mean curvature motion is seen to be the application of the diffusion process to the geometry itself, rather than the diffusion of function values defined over the geometry.
One of the shortcomings of using diffusion filtering in its basic form (7.2) for surface texture purposes is that it has a tendency to smooth out sharp edges as well as the general noise. Perona & Malik (1987) proposed a method to overcome this by introducing nonlinearity into equation (5.1) in the form of a diffusivity function. The diffusivity function is usually chosen to yield very small values for regions that have large differences between local surface height values. The diffusivity function acts like a weight function for the diffusion process, with the result that diffusion is reduced in the vicinity of edges, while general noisy regions are smoothed normally.
Another development is the application of anisotropic diffusion. This allows the diffusion to be direction dependent via the introduction of a diffusivity tensor (Weickert 1998).
8. Discrete differential geometry
In general, for most practical applications, the geometry we are looking to diffuse takes the form of a three-dimensional polygon mesh formed from discrete data points. In such cases, there are a number numerical methods available for approximating both the Laplace–Beltrami operator and solutions of the diffusion equation at discrete time values. This is particularly true when the polygon mesh is composed of triangles.
(a) Discrete Laplace–Beltrami operators
Most discrete approximations to the Laplace–Beltrami operator involve triangular meshes, although Liu et al. (2008) describes a method for use with quadrilateral meshes. We confine our interest to the triangular mesh approximations proposed by Meyer et al. (2002) and by Desbrun et al. (1999), mainly owing to the study of Xu (2004), who establishes the convergence of these schemes when the surface triangulation in parameter space satisfies certain criteria (described later). Some notation required for the description of discrete Laplace–Beltrami approximations is now introduced. This notation is consistent with that of Xu and Meyer.
A triangular mesh M is considered to be composed of N measured data points that provide a discrete approximation of a continuous but unknown underlying surface . Edges are defined as straight line segments joining two neighbouring vertices pi and pj and are denoted by [pi,pj]. Triangles are denoted by [pi,pj,pk], where pi,pj,pk are the vertices of the triangle. N(i) is the set of all points pj that are vertices of triangles containing pi (members of N(i) are called one-ring neighbours of pi). The cardinality of the set N(i) is referred to as the valence of vertex pi. Using this notation, the Desbrun and Meyer approximations for the Laplace–Beltrami operator applied to a function are, respectively, given by8.1and8.2Here, A(pi) is the area of all triangles that contain vertex pi, and the angles αij and βij are the two angles opposing the edge pi,pj of the two triangles that share edge pi,pj. AM(pi) is the area of the region formed by joining the circumcentres of all triangles containing vertex pi. In cases where a triangle is obtuse, the midpoint of the edge opposite the obtuse angle is used instead of the circumcentre (Meyer et al. 2002). Figure 2 provides a visual representation of the angles αij and βij.
(b) Criteria for convergence
Xu (2004) has shown that the Meyer (8.2) and Desbrun (8.1) approximations are shown to converge to the Laplace–Beltrami operator at a point, subject to certain conditions on the positioning of each point and its one-ring neighbours in parameter space. The convergence is defined in the sense that the theoretical Laplace–Beltrami operator is the limiting value of these approximations calculated on increasingly fine subdivisions of the initial triangular mesh. The first condition for convergence is that the vertex for which we are calculating the approximation has a valence of 6. The second condition is that in parameter space, each vertex and its one-ring neighbours form a sheared hexagon. Formally, if the calculation point has six one-ring neighbours , each with corresponding parameter values qc,q1,…,q6 where , such that pc=x(qc), , then convergence will occur if8.3where index arithmetic is performed modulo 6. The proofs of convergence hold for the mean curvature H(p) values as well as for the Laplace–Beltrami operator of a function . Further details and proofs of the convergence for these approximations can be found in Xu (2004).
(c) Numerical solutions of the diffusion equation
To solve the diffusion equation at a point in time, we have used the semi-implicit Euler scheme as described by Desbrun et al. (1999). Assuming that the surface is evolving in time, we denote the matrix of surface points at time t by . Row i of this matrix contains the space coordinates of vertex pi. For a time step of δt from an initial time t0, we approximate the time derivative using a finite difference8.4Now the location of the surface points as a result of the diffusion process (5.1) has an approximate solution P(t+δt) at time t+δt given by the linear system8.5The matrix represents the discrete approximation of the Laplace–Beltrami operator using either equation (8.2) or (8.1). For the Meyer et al. approximation (8.2), Δ has element i,j given by8.6For the Desbrun et al. (8.1) approximation, Δ has element i,j given by8.7From these matrices, it can be seen that the discrete Laplace–Beltrami operator for vertex i at time t is given by the ith row of the matrix product ΔP(t). The equivalent formulation can be used to solve the diffusion equation for functions f on the surface . In such cases, the matrices are replaced with , which has the function value f(pi) at time t as its ith element. In either case, the system can be solved using standard techniques for linear systems, although given the sparse nature of the matrix Δ, methods for solving sparse linear systems may be more appropriate. This is of particular relevance for the analysis of a large number of data points. Desbrun et al. (1999) propose the use of a preconditioned bi-conjugate gradient descent methods to provide an iterative solution to the spare linear system.
9. Proposed method of application
Application of mean curvature motion as described above using both Meyer and Desbrun forms of the discrete Laplace–Beltrami operator was found to be unsatisfactory for the purposes of surface metrology. Results obtained from these methods using simulated and experimental data in the form of surface meshes with boundary were found to have an unacceptable amount of shrinkage and drift in space. Owing to the nature of surface metrology, precision is of major importance and as a consequence the unpredictable results of mean curvature motion were considered to be unsuitable for data analysis.
As it is an important requirement for the filtered data to remain close to the initial set of data points, an alternative approach is proposed. Instead of performing diffusion of the geometry itself, a least-squares approximation is fitted to the surface data to provide a reference surface, and diffusion filtering of the approximation residuals is performed on the geometry of the reference surface. Firstly, a least-squares NURBS surface approximation is fitted to the initial set of N measurement points . This reference surface is denoted by (as defined in equation (3.1)). Points (where ) on the surface are then calculated for each measurement point ri. The least-squares approximation residual ei at this point is also calculated. In this way, samples of a function defined by f(pi)=ei are generated. The idea is to then perform diffusion filtering of the residuals ei on the new mesh obtained from the reference points pi on the reference surface using the methods described in §8. Triangulation of the points qi that form an approximate rectangular grid in parameter space is performed in the same way as proposed for mean curvature motion: form two triangles from each quadrilateral in the mesh by joining its lower left and top right vertices. This again ensures a valence of six for all interior vertices, and very closely approximates the convergence condition (8.3). By choosing to diffuse the least-squares residuals on the reference surface, the data are unable to drift too far from its initial configuration. Additionally, it is reasonable to assume (in the case of a good surface approximation) from statistical theory that the least-squares residuals ei will be approximately normally distributed with mean zero (ei∼N(0,σi)). In such a case, it would be expected that over time, the diffusion process would result in the residuals converging to their mean value . Assuming the validity of this statistical reasoning and provided that surface is a good approximation to the initial data, it is reasonable to expect to be very close to zero. As a consequence, the limiting surface of this diffusion process will be very close to the reference surface .
(a) Surfaces without boundaries
Most of the literature on diffusion filtering deals with the case of Euclidean geometry or more general geometries where the data in question are perceived as measurements from an underlying manifold that has no boundary. In cases where the manifold has no boundary, the triangulation can be closed, and so there are no issues with evaluation of curvature at each point as every vertex has a complete set of one-ring neighbours. The use of equations (8.1) and (8.2) for calculation of curvature is not particularly satisfactory at edge vertices, and experimentation with these methods has not yielded good results at boundaries in such cases. At such a time that surface measurement data take the form of a closed mesh without boundary, this will not be an issue. Currently, most surface measurement data take the form of surface sheets with a clear boundary and so a method of dealing with boundary vertices is required for a satisfactory data analysis.
(b) Surfaces with boundaries
Because of the undesirable results that we obtain for the Laplace–Beltrami operator at edge vertices, they are used only in the calculations (8.2) and (8.1) for the Laplace–Beltrami operator at interior vertices. They are included in the calculations, but their actual diffused residual values are disregarded for analysis purposes. This ensures that only vertices that have a valence of six are actually filtered to ensure the reliability and accuracy of the filtered function values. Providing that the measured surface is initially sampled at a sufficiently large amount of points, this approach should not result in the loss of a large percentage of data. Xu (2004) suggests that for diffusion of surfaces with a boundary, a possible way of dealing with the problem is to add extra surface points around the boundary to provide an outer layer of triangles. This results in the initial boundary vertices having an artificial valence of six. For surface fairing, where high accuracy may not be the chief consideration, this may be an acceptable approach, but in a field such as surface metrology where data accuracy and reliability is of high importance, this seems to be an unreliable way to proceed. For the method of filtering residuals on a reference surface, this problem would be exacerbated by the need to simulate a residual value for each new artificial surface point. This approach provides no means of assessing how much inaccuracy is introduced by the simultaneous extrapolation of geometry and simulation of function values. For these reasons, at this point in time, the exclusion of boundary vertices for filtering purposes is proposed.
10. Computing simulation studies
For the purposes of testing these methods, a simple NURBS surface (degree (4,4) with 6×6 control points) was created as the reference surface and is shown in figure 3. For the surface points, a uniform rectangular mesh was created in the parameter domain, and this was then triangulated in the manner described in §9. In order to simulate a realistic surface roughness sample, the reference surface residuals were selected as a subset of a stochastic model that generated a roughness surface (derived from the internet-based surface metrology algorithm testing system of the National Institute of Standards and Technology; NIST 2010). These residuals were then filtered using the diffusion equation as described previously using the Meyer form of the discrete Laplace–Beltrami operator (8.2). Figures 4–7 illustrate the results of the diffusion process at increasing values of the time parameter t that correspond to common choices of the cut-off wavelength λc. Owing to the fact that the surface form change is much larger than the residuals, the surface residuals are plotted against their values in parameter space in order to enhance the visibility of the filtration effects. As mentioned in §9, the boundary residuals are not diffused, and so only the internal vertex residuals are illustrated in figures 4–7.
11. Experimental studies
After the initial application of these methods to simulated data, diffusion filtering was performed on real surface measurement data. The data were obtained from coordinate measuring machine (CMM) measurement of a portion of a worn knee replacement component. The data are composed of 60×35 points taken on an approximately rectangular grid. The approximating reference surface is a cubic NURBS surface with 14×17 control points and is illustrated in figure 8. As with the simulated data, the diffusion of the fitted surface residuals was implemented using the Meyer et al. approximation of the Laplace–Beltrami operator (8.2). Figures 9–12 show the results of diffusion of the residuals at increasing values of the time parameter. As with the simulated data results, the residuals are plotted against the parameter values of the surface reference points so that the effects of filtration are clearly visible. Boundary vertices are again omitted from these figures for the same reasons as mentioned in §10.
The results show that a potential way forward for filtration of freeform surface metrology data is to use the established field of PDE filtering. In particular, diffusion filtering can be applied to obtain a Gaussian scale space for freeform surfaces owing to its theoretical link with the linear Gaussian filter as described in §5. As PDE filtering is a wide and active area of research, the theory of these methods is well discussed in other research fields, such as image analysis. This paper is a pilot study and is intended to provide a useful starting point for solving freeform filtration difficulties by using PDEs. The proposed use of a reference surface to provide a manifold on which to apply the diffusion equation provides a method of ensuring that the filtered data do not deviate excessively from the initial data. This is an important requirement for analysing surface measurement data, and overcomes the undesirable mesh area shrinkage and spatial drift effects that were exhibited by experimentation with mean curvature motion. It is hoped that further research into these methods will yield a way to apply mean curvature motion in such a way that these undesirable aspects can be removed. Much more research into this area is required, but the results of the initial experimentation with these methods is encouraging.
13. Future work
A key area of future research will be the application of this methodology to a large number of different types of freeform surface data, and to establish whether there are any particular surfaces for which the method yields very good or poor results. Other areas of interest will be development of ways of speeding up the algorithms, and an investigation into the propagation of uncertainty through the filter.
Other definitions for the discrete Laplace–Beltrami operator are available, and a deeper investigation of these would seem appropriate, particularly mesh-free definitions. Experimentation with different types of (and higher order) geometric flow is another area that is likely to be worthy of attention. So far, this paper has only considered linear isotropic diffusion as a means of filtering the data, but experimentation with anisotropic diffusion or other edge-enhancing types of diffusion may lead to superior results, particularly in the analysis of structured surfaces. Finally, it would be desirable to try and establish a theoretical link between the cut-off wavelengths of the linear areal Gaussian filter and the time step in the evolution of the diffusion equation.
X.J. gratefully acknowledges the Royal Society under a Wolfson-Royal Society Research Merit Award, and the authors gratefully acknowledge the European Research Council for its ‘Ideal Specific programme’ ERC-2008-AdG 228117-Surfund.
- Received June 13, 2010.
- Accepted August 13, 2010.
- This journal is © 2010 The Royal Society