Robust surface profile envelope estimation using a beam-surface contact model

Envelope lines of surface profiles have found application in the function oriented analysis of technical surfaces [1–3]. For example, plateau-like surfaces where the characteristics of the valleys are critical for the technical functionality [2]. These envelope lines are commonly constructed using morphological filters [4] or motif methods [5]. While none of these methods use contact mechanic models, they have been successfully used to describe the contact behavior of surfaces [1]. However, they do suffer from some drawbacks. The morphological filters are sensitive to isolated peaks and produce boundary artifacts [6, 7]. The adequate selection of the structuring element for the filtering of textured surfaces is also the topic of recent research [8]. The motif methods are sensitive to lateral shifts and the sampling direction [5].

In order to provide a robust starting point for the function oriented analysis of technical, plateau-like contact surfaces, the estimated envelope should possess several characteristics. Isolated narrow peaks should be largely ignored, since they likely will not provide sufficient bearing for a contact partner and will quickly wear down or are simply particles of dirt. Valleys in the surface should, depending on the desired scale of the extracted envelope line and valley width, either be also ignored or incorporated by the envelope line [7]. It is also desirable for the envelope line to be smooth.

In this paper we describe a method that estimates such an envelope line based on the bending line of a beam pressed against the surface profile. The notion of a beam in the context of surface profile filtration is closely associated with spline filters [9]. In [9], the mathematical equations for the spline filter and its properties were investigated first, and then the mechanical model of the beam deformed by springs (see figure 1(a)) was presented to provide some intuition about the spline filter.

Zoom In Zoom Out Reset image size

We approached the development of the proposed method differently. We wanted a method to calculate envelope lines that are as smooth, robust against outliers and flexible in their shape as the mean lines of the spline filter and its robust variant. Then we derived the necessary changes in the mechanical model of the spline filter that would change the resulting beam bending line from a mean line to an envelope line. The resulting mechanical model is drawn in figure 1(b), and the corresponding equations as well as their properties are analyzed in section 2.

Looking at the mechanical models of the classical smoothing spline in figure 1(a) and the proposed mechanical model in 1(b), the differences between these two methods become apparent. In both cases, the interaction between the discretized profile points and the beam is modeled as springs. But the springs of the spline filter model are attached to the beam and the profile points. The spring forces either push the beam up or down, depending on the height of the profile point relative to the deformed beam, and the resulting beam bending line will consequently approximate the local mean of the profile heights.

In the proposed beam-surface method the springs are attached neither to the profile points nor to the beam. Instead, the lower spring end is attached to a virtual ground below the profile and the upper end is unattached. The length of each spring in the uncompressed state is such, that it reaches exactly the corresponding profile point from the virtual ground. This is the case for the outermost left and the two outermost right springs. If a profile point lies above the deformed beam, its corresponding spring is in contact with the beam and pushes the beam upwards. The only downwards force is provided by the line load. Therefore, a profile point does not affect the beam as long as it does not penetrate the profile at the given point.

In this section we will derive the mathematical equations corresponding to the mechanical model illustrated in figure 1(b), how to solve the resulting system of equations and an intuitive parameterization of said equation system.

2.1. Equations

The model consists of an Euler-Bernoulli beam with flexural rigidity EI and a perpendicular line load q. In this paper, the load is assumed to be distributed uniformly along the beam but other distributions are usable as well.

The beam is discretized using two knot elements, with three degrees of freedom per knot. These degrees of freedom correspond to the beam bending line w(x) position z_i , slope φ_i and curvature χ_i at the knot positions. The six degrees of freedom of a single element u_1,2 are then

$$\begin{eqnarray}{{\bf{u}}}_{\mathrm{1,2}}={\left[\begin{array}{cccccc}{z}_{1} & {\varphi }_{1} & {\chi }_{1} & {z}_{2} & {\varphi }_{2} & {\chi }_{2}\end{array}\right]}^{T}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&w({x}_{i})={z}_{i},\,\displaystyle \frac{{\rm{d}}w({x}_{i})}{{\rm{d}}x}={\varphi }_{i},\,\displaystyle \frac{{{\rm{d}}}^{2}w({x}_{i})}{{\rm{d}}{x}^{2}}={\chi }_{i}.\end{eqnarray} \tag{ 2 }$$

According to [10], the linearized stiffness matrix K_1,2 for such an element is

$$\begin{eqnarray}{{\bf{K}}}_{\mathrm{1,2}}=\displaystyle \frac{{EI}}{35{l}^{3}}\left[\begin{array}{cccccc}600 & 300l & 15{l}^{2} & -600 & 300l & -15{l}^{2}\\ 300l & {192}^{2} & 11{l}^{3} & -300l & 108{l}^{2} & -4{l}^{3}\\ 15{l}^{2} & 11{l}^{3} & 3{l}^{4} & -15{l}^{2} & 4{l}^{3} & \displaystyle \frac{1}{2}{l}^{4}\\ -600 & -300l & -15{l}^{2} & 600 & -300l & 15{l}^{2}\\ 300l & 108{l}^{2} & 4{l}^{3} & -300l & 192{l}^{2} & -11{l}^{3}\\ -15{l}^{2} & -4{l}^{3} & \displaystyle \frac{1}{2}{l}^{4} & 15{l}^{2} & -11{l}^{3} & 3{l}^{4}\end{array}\right],\end{eqnarray} \tag{ 3 }$$

where EI is the flexural rigidity and l the length of the beam element, and a constant line load q is distributed among the coordinates as

$$\begin{eqnarray}{{\bf{q}}}_{\mathrm{1,2}}=\displaystyle \frac{q}{120}{\left[\begin{array}{cccccc}60 & 12l & {l}^{2} & 60 & -12l & {l}^{2}\end{array}\right]}^{T}.\end{eqnarray} \tag{ 4 }$$

The static equation of the mechanical model depicted in figure 2 can then be expressed as

$$\begin{eqnarray}&&{{\bf{K}}}_{\mathrm{1,2}}{{\bf{u}}}_{\mathrm{1,2}}={{\bf{q}}}_{\mathrm{1,2}}+\left[\begin{array}{c}\sigma ({v}_{1}-{z}_{1})({v}_{1}-{z}_{1})c\\ 0\\ 0\\ \sigma ({v}_{2}-{z}_{2})({v}_{2}-{z}_{2})c\\ 0\\ 0\end{array}\right]\end{eqnarray} \tag{ 5 }$$

where σ is the step function, v_i is the z position of the profile and c the stiffness of the springs. These equations can easily be extended to a profile of N points, with a knot at each sample point. The left-hand side of the resulting equation system is the sum of the beam stiffness K, a banded symmetric positive definite (SPD) matrix, and a diagonal positive semi-definite matrix C representing the stiffness of the bearing. The sum will therefore remain a banded symmetric positive definite matrix. The right-hand side is a sum of the line load q and a term representing the initial position of all material points in contact c.

$$\begin{eqnarray}&&\left({\bf{K}}+{\bf{C}}({\bf{u}})\right){\bf{u}}={\bf{q}}+{\bf{c}}({\bf{u}})\end{eqnarray} \tag{ 6 }$$

In between the discontinuities at z_i = v_i , where a knot makes or breaks contact with the material, the static equation can be expressed as a classic linear system of equations Au = b. Each making or breaking of contact can be represented as a rank-1 update or downdate of A. Therefore, we cannot directly solve the equation system since the local A, corresponding to the current solution estimate, will be related to the final A by an up- or downdate of unknown rank.

Zoom In Zoom Out Reset image size

2.2. Solver

An intuitive approach to solve the equation system (6) is an iterative solver and since A is guaranteed to be SPD, the conjugate gradient method is a promising start [11]. During each iteration a step size β_k and update direction p_k is calculated such that

$$\begin{eqnarray}&&{{\bf{u}}}_{{\bf{k}}+1}={{\bf{u}}}_{{\bf{k}}}+{\beta }_{k}{{\bf{p}}}_{{\bf{k}}}\end{eqnarray} \tag{ 7 }$$

conforms to

$$\begin{eqnarray}&&{\parallel {\bf{b}}({{\bf{u}}}_{{\bf{k}}+1})-{\bf{A}}({{\bf{u}}}_{{\bf{k}}+1}){{\bf{u}}}_{{\bf{k}}+1}\parallel }_{2}\lt {\parallel {\bf{b}}({{\bf{u}}}_{{\bf{k}}})-{\bf{A}}({{\bf{u}}}_{{\bf{k}}}){{\bf{u}}}_{{\bf{k}}}\parallel }_{2}.\end{eqnarray} \tag{ 8 }$$

The proposed algorithm uses the preconditioned conjugate gradient method to calculate the update direction p_k and initial step size guess α_k . Since there is no guarantee that the solution u_k + α_k p_k satisfies the monotonicity condition equation (8), we also employ a backtracking line-search algorithm with the Armijo-Goldstein condition [12] to compute the conforming step size β_k (Algorithm 1). Starting from the maximum step size α_k , we will check if the ratio of the true and the linearized improvement is higher than a control parameter c ∈ (0, 1]. If not, we multiply the step size with the second control parameter τ ∈ (0, 1) and repeat until linearized and true improvement are close enough. Since the linearization is exact up to round-off errors between the discontinuities and the problem is convex, there is guaranteed to be a step size that reduces the residual.

Algorithm 1. Backtracking line search.

Input: ${{\bf{u}}}_{{\bf{k}}},{{\bf{p}}}_{{\bf{k}}},{\alpha }_{k},\tau ,c$

Output: ${\beta }_{{\bf{k}}}$

begin

j = 0

${r}_{0}={\parallel {\bf{b}}({{\bf{u}}}_{{\bf{k}}})-{\bf{A}}({{\bf{u}}}_{{\bf{k}}}){{\bf{u}}}_{{\bf{k}}}\parallel }_{2}$

repeat

${\beta }_{k}={\tau }^{j}{\alpha }_{k}$

${{\bf{u}}}_{{\bf{k}}+1}={{\bf{u}}}_{{\bf{k}}}+{\beta }_{k}{{\bf{p}}}_{{\bf{k}}}$

$\hat{r}={{\bf{r}}}_{0}-{\parallel {\bf{b}}({{\bf{u}}}_{{\bf{k}}})-{\bf{A}}({{\bf{u}}}_{{\bf{k}}}){{\bf{u}}}_{{\bf{k}}\,+\,1}\parallel }_{2}$

$r={{\bf{r}}}_{0}-{\parallel {\bf{b}}({{\bf{u}}}_{{\bf{k}}+1})-{\bf{A}}({{\bf{u}}}_{{\bf{k}}+1}){{\bf{u}}}_{{\bf{k}}\,+\,1}\parallel }_{2}$

$j=j+1$

until $r/\hat{r}\geqslant c$

return ${\beta }_{k}$

end

For good convergence behavior, we precondition our local equation system with an incomplete Cholesky decomposition [11]. This is not a very expensive operation since for banded SPD matrices of bandwidth m and size n × n the full Cholesky decomposition has a complexity of ${ \mathcal O }({m}^{2}n)$.

Algorithm 2. Beam solver.

Input: ${\bf{A}}({\bf{u}}),{\bf{b}}({\bf{u}}),{{\bf{u}}}_{0},{r}_{\max }$

output: ${{\bf{u}}}_{{\bf{k}}}$

begin

k = 0

${{\bf{u}}}_{{\bf{k}}}={{\bf{u}}}_{0}$

${{\bf{r}}}_{{\bf{k}}}={\bf{b}}({{\bf{u}}}_{{\bf{k}}})-{\bf{A}}({{\bf{u}}}_{{\bf{k}}}){{\bf{u}}}_{{\bf{k}}}$

${r}_{{\bf{k}}}={\parallel {{\bf{r}}}_{{\bf{k}}}\parallel }_{2}$

while ${r}_{k}\gt {r}_{\max }$ do

/∗ incomplete Cholesky decomposition ∗/

${{\bf{L}}}_{{\bf{k}}}=\mathrm{ichol}\,\left({\bf{A}}({{\bf{u}}}_{{\bf{k}}})\right)$

${{\bf{L}}}_{{\bf{k}}}{{\bf{y}}}_{{\bf{k}}}={{\bf{r}}}_{{\bf{k}}}$ /∗ forward substitution ∗/

${{\bf{L}}}_{{\bf{k}}}^{T}{{\bf{p}}}_{{\bf{k}}}={{\bf{y}}}_{{\bf{k}}}$ /∗ back substitution ∗/

${\alpha }_{{\bf{k}}}=\tfrac{{{\bf{r}}}_{{\bf{k}}}^{T}{{\bf{p}}}_{{\bf{k}}}}{{{\bf{p}}}_{{\bf{k}}}^{T}{{\bf{Ap}}}_{{\bf{k}}}}$

${\beta }_{{\bf{k}}}=\mathrm{linesearch}\,\left({{\bf{u}}}_{{\bf{k}}},{{\bf{p}}}_{{\bf{k}}},{\alpha }_{{\bf{k}}}\right)$

${{\bf{u}}}_{{\bf{k}}+1}={{\bf{u}}}_{{\bf{k}}}+{\beta }_{k}{{\bf{p}}}_{{\bf{k}}}$

${{\bf{r}}}_{{\bf{k}}+1}={\bf{b}}({{\bf{u}}}_{{\bf{k}}+1})-{\bf{A}}({{\bf{u}}}_{{\bf{k}}+1}){{\bf{u}}}_{{\bf{k}}\,+\,1}$

${r}_{k+1}={\parallel {{\bf{r}}}_{{\bf{k}}+1}\parallel }_{2}$

$k=k+1$

end

return ${{\bf{u}}}_{{\bf{k}}}$

end

2.3. Normalization and parameterization

Two issues need to be addressed with this approach. First, our mechanical model is based on the linearized beam equations. Therefore, the results will mismatch the intuitive expectations in case of large slopes in the fitted beam. This can be prevented by normalizing the input profile data. Second, it is not obvious how to choose the parameters EI, c, and q to achieve the desired filter effect. We therefore introduce two new parameters, the nesting index NI and the average penetration depth z_p , from which EI, c, and q can be calculated.

One of the underlying assumption of the linearized beam equation is that the slope of the beam is small enough for the following simplification to be valid

$$\begin{eqnarray}&&\displaystyle \frac{{\chi }_{i}}{\sqrt{{\left(1+{\varphi }_{i}^{2}\right)}^{3}}}\approx {\chi }_{i},\quad \mathrm{if}\ {\varphi }_{i}^{2}\ll 1.\end{eqnarray} \tag{ 9 }$$

Verifying this requires knowledge of the solution to the beam equation which is not available at the beginning of the algorithm, when the normalization is performed. We therefore calculate the scaling factor such that the squared median slope of the input profile data equals 0.01. The average squared slope of the beam will be smaller than the average squared slope of the input profile data and therefore equation (9) will be a valid approximation for most profiles. The validity of the approximation can be verified using the final beam bending line and if necessary a different scaling factor can be selected.

As stated above we introduce the nesting index NI and the average penetration depth z_p to substitute the flexural rigidity EI, the beam line load q and the bearing stiffness c. Since we substitute three free parameters with two free parameters one of the three parameters needs to be fixed to get an unambiguous relation. We choose a fixed value of 1 × 10⁻⁴N mm² for EI and continue to derive the mapping of z_p and NI to the ratios EI/q and q/c. Numerical simulations have shown that this value for EI mostly leads to sufficiently well conditioned equation systems.

The average penetration depth results from the static equivalence of the beam line load and the elastic stress of the springs representing the profile points. We define z_p as the average difference between the input profile and the resulting beam bending line for the case of the input profile being of constant height. The ratio between the line load q and the spring bearing stiffness c is therefore computed as

$$\begin{eqnarray}&&\displaystyle \frac{q}{c}=\displaystyle \frac{-N}{(N-1)l}{z}_{p}.\end{eqnarray} \tag{ 10 }$$

The flexural rigidity EI of the beam acts like a low-pass filter, increasing EI leads to smoother bending lines. An increase in the line load q in turn results in a bending line that follows the input profile closer. In accordance to [13], we define the nesting index NI as the wavelength of a pure cosine input profile whose beam bending line amplitude is 6dB smaller than the input profile amplitude. The amplitude of the input cosine is normalized according to the previously described scheme. The nonlinear relationship between NI and EI/q for several penetration depths z_p is depicted in figure 3.

Zoom In Zoom Out Reset image size

Some beam bending lines for different NI and z_p are shown in figure 4. As expected, the smoothness of the beam bending line rises with increasing NI. For small z_p the bending line rests on the peaks of the profile. With increasing penetration depth, more points of the profile come in contact and the influence of outliers in positive z direction is reduced.

Zoom In Zoom Out Reset image size

In this section, we will present and discuss some results from applying the algorithm to different profiles from sintered surfaces. These surfaces are typically plateau-like with deep notches. We will also compare it with commonly used morphological operators and present another possible application of the beam-surface contact model.

Four examples of sintered surface profiles with beam envelopes are shown in figure 5. The beam bending line is relatively unaffected by the notches but follows trends in the profile that are longer than the nesting index of NI = 3 mm. There are several isolated peaks in the profiles that do not visibly affect the envelopes.

Zoom In Zoom Out Reset image size

Commonly, morphological operators are used to extract the envelope from a profile [1, 4]. The closing operator, which is a dilation followed by an erosion, is usually employed for this task. The closing operation is very sensitive to outliers in positive z direction, as demonstrated in figure 6. To reduce the sensitivity against these outliers, the closed profile can subsequently be filtered using an opening operator with a larger structuring element [6].

Zoom In Zoom Out Reset image size

Another interesting result is obtained by sweeping over different z_p and plotting the percentage of bearing profile points over z_p . This simulates the contact behavior of the surface depending on the contact pressure, since z_p is proportional to the applied line load q. As illustrated in figure 7, the results look similar to the Abbott curve used in the characterization of surface profiles. At small penetration depths of z_p = 10 × 10⁻⁴mm, a much larger percentage of the profile area of Profile 2 and Profile 3 is bearing the load compared to Profile 1 and Profile 4. At deeper penetration depths the difference in the bearing contact area decreases again. Unlike the original Abbott curve, which is a cumulative histogram of profile z-values, the neighborhood of each profile point is also taken into account in this approach.

Zoom In Zoom Out Reset image size

This paper describes a new method to construct an envelope for a discretized surface profile based on a contact mechanic model between a beam and a rough surface. Such an envelope is often used as a starting point for the functional analysis of the profile like assessing the tribological behavior of the surface or detecting surface defects.

The method can be tuned using the two parameters z_p and NI. The average penetration depth z_p determines the ratio between the line load pressing the beam into the surface and the stiffness of the profile. The nesting index NI determines the smoothness of the envelope and is closely related to the cutoff wavelength of a linear filter. Their influence on the resulting envelope has been shown in figure 4. Due to the intuitive nature of the parameters it should be straightforward to determine suitable parameterization for different profile types.

The behavior of the proposed method compared to morphological methods is exemplified in figure 6. We will highlight three aspects: Robustness against outliers, end effects and envelope shape.

The sensitivity of the closing operation to isolated outliers and its reduction due to a successive opening operation is apparent. The beam bending line does not show any comparable deflections due to narrow peaks. This robustness could be even further increased by incorporating an upper limit for the spring force at each point, representing a plastic deformation of the profile in our mechanical model. The shape of the used structuring element is still visible in the envelope of the morphological filter. Since the underlying function shape of the proposed method is a cubic spline, the resulting envelope is flexible in its shape while remaining smooth. Zooming in on the right end of the plot in figure 6, we notice a sharp drop in the envelopes generated by the morphological methods. This is due to the profile ending in a valley and because the recommended boundary conditions from [4] have been applied. The natural boundary conditions used by the proposed method do not cause such extreme end effects and the entire measurement length can be used for evaluation.

Surface characterization based on envelope-lines are seen as an important complement to surface characterizations based on mean-line filters. We introduce a new method for envelope estimation based on a beam-surface contact model that shares some benefits of robust mean-line filters like the robust spline filter or the robust Gaussian regression filter [14]. These benefits are robustness against outliers, small end distortions and a high flexibility regarding the shape. In a future paper the relationship of the proposed method to the robust spline filter will be investigated in more depth.

Another possible application of the proposed method next to constructing envelopes was hinted at in figure 7. By continuously increasing the average penetration depth and plotting the corresponding percentage of bearing profile points we receive a curve similar to the Abbott or bearing area curve. The main difference being that the curve shown in figure 7 takes the neighborhood of each point into account while the Abbott curve solely uses the height of each profile point individually. Since the proposed method is based on a contact mechanic model, it might provide the basis for a more accurate assessment of the surface bearing properties compared to the classical Abbott curve. We will further explore this idea in a subsequent paper.

All data that support the findings of this study are included within the article (and any supplementary files).