Available online at www.sciencedirect.com

ScienceDirect

Procedia CIRP 43 (2016) 29 - 34

14th CIRP Conference on Computer Aided Tolerancing (CAT)

A sampling strategy based on B-wavelets decomposition

Luca Pagania*, Paul J. Scott

aEPSRC Centre for Innovative Manufacturing in Advance Metrology, Centre for Precision Technologies (CPT), School of Computing and Engineering, University

of Huddersfield, Huddersfield, HD13DH, UK * Corresponding author. Tel.: +44-1484-471284; fax: +44-1484-472161. E-mail address: l.pagani@hud.ac.uk

Abstract

Finding optimal sampling strategy is a central problem in surface reconstruction. In this paper a sampling strategy on different scales of a free-form surface is proposed. Non-Uniform Rational B-spline (NURBS) surfaces are used to represent the nominal shape; the decomposition is performed through the orthogonal B-wavelet basis. A sampling based on each level of the shape decomposition is then proposed to detect the points of interest that reflect the variations between different levels. The number of the samples is selected according to the complexity of the shape, which is represented by the number of the internal knots.

© 2016 The Authors.PublishedbyElsevier B.V. Thisis an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

Peer-review under responsibility of the organizing committee of the 14th CIRP Conference on Computer Aided Tolerancing Keywords: Sampling; B-splines; B-wavelets; Decomposition; Curve Reconstruction; Surface Reconstruction.

1. Introduction

Finding optimal sampling strategy is a central problem in surface reconstruction as well as in surface inspection. In this paper we propose a sampling strategy based on the decomposition of the curve or surface in different scales. Since NonUniform Rational B-spline (NURBS) curves and surfaces are widely used in Computer Aided Geometric Design (CAGD), we propose to directly use these functions to represent the nominal shape. We firstly propose to decompose the nominal object at different approximation scales with the orthogonal B-wavelet decomposition [1,2], a sampling according to the differences between each level can then be computed. The NURBS functions are used instead the triangular mesh in order to maintain the exact geometry of the object. Representing the real object through its NURBS representation has several advantages also in finite element method simulation [3]. One of the advantages of using the B-wavelets decomposition instead of using the discrete wavelet transform is that the number of the levels of the B-wavelets functions depend on the number of the internal knots of the associated B-splines function, while those of the discrete transform depend on the number of sampled points, and this number can increase rapidly if the sampling rate is high. On the contrary the reconstruction with the B-splines function can bound the number of decomposition levels. A profile sampling based on the decomposition with the discrete wavelets transform has been proposed in [4]. In this paper the discrete wavelet transform is used to obtain the curvature information of a sig-

nal. The samples are then allocated according to the curvature distribution in different regions. In the proposed method we use all the information coming from the different approximation scales.

The paper is structured as follows: in Section 2 we introduce the NURBS functions and the B-wavelets transformation, in Section 3 we present the decomposition algorithm, in Section 4 we explain the sampling method and in Section 5 the reconstruction performance.

2. B-wavelets decomposition

Before introducing the B-wavelets decomposition [1], we briefly recall the construction of B-splines and NURBS functions [5]. Given a non decreasing sequence of knots t = (ii)n=Î+\ the B-splines base BUr(x) = BUrt(x) of degree r ^ 0 with support [ti, ti+1] is defined by the De Boor's algorithm by

Bi,r(t) = t Bi,r-1(t) + -t+r—LBi+1,r-1(t) (1)

ti+r-1 ti

■ ti+1

Bi,0(t)

1 if ti < t < ti+1 0 otherwise.

2212-8271 © 2016 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.Org/licenses/by-nc-nd/4.0/).

Peer-review under responsibility of the organizing committee of the 14th CIRP Conference on Computer Aided Tolerancing doi:10.1016/j.procir.2016.01.028

A r-th degree B-splines curve is defined as

r(t) = ^ Bi ,r№ Pi tl ^ t ^ trt+r+l

w(t) = Yj Wi Bi,d(t).

where pt e Rd are the control points.

In order to approximate conics, the r-th degree NURBS curve has been introduced, it is defined as

r(t) ■

Z?=o WiBir (t) Pi 2^=0 WiBi,r(t)

tl < t < tn+d+1

where wi are the weights of the i-th B-splines basis function.

A NURBS surface is defined through the tensor product of two NURBS basis functions in two orthogonal directions u and v. A NURBS surface of degree p in u direction and q in v direction is a piecewise rational function

r(u, v)

£?=oZ7=0 WijBip(u) Bjqq(v) Pij

= 2?=o j wijBip(u) Bjq(v) (5)

with ui < u < u?+p+i, vi < v < v?+p+i

The space Vd can therefore be written as the direct sum Vd = Vd-1 ffi Wd-1 where Wd-1 is the orthogonal component of Vd-1 in Vd .By iterating the decomposition d time we get

Vd = Wd-1 ® Wd-2 ® ••• ® Wo ® Vo.

The functions in the spaces W, are called wavelets, and the corresponding spaces are called wavelets spaces. B-wavelets are wavelets function with minimal support [1]. In the remaining of the paper we will refer to the B-wavelets basis function with the letter i^(t).

We now focus on the decomposition of the space V into the spaces Vi-1 and Wi-1. The space V is spanned by the B-spline basis functions identified by the knots vector while the knots vector of the B-splines in the space Vi-1 is ti-1 that is build by deleting one or more internal knots from the vector t,.

Let = (Br>1,..., Br,ni)T denote the vector of the B-spline basis functions of the space V, and = (i^,..., fni-nil)T the vector of the B-wavelets basis function of the space Wi, a points of the B-splines curve of the space Vi can be expressed as

where {ptj} are the control points, wtj are the weights of the product of the B-splines basis B;(w) and Bj(v), u and v are the knots vectors of the B-spline basis in u and v direction.

The non-uniform B-wavelet transform has been introduced in [1] and it has been extended to the NURBS space in [2]. The B-wavelet transform is a mathematical tool designed to decompose a B-spline function belonging to the space Vd through a projection in two spaces Vd-1 and Wd-1. The space Wd-1 consist of functions in Vd that are orthogonal to the functions of the space Vd-1

ri(u) = </>] (u) Ci

where Ci = [piA,..., Pinif (Ci = [wn pa ..., Win Pi,n]T) is the vector of the B-spline (NURBS) curve control points of the i-th level. It is possible to rewrite the functions in the vectors as a linear combinations of the functions and tfri_1 as [1]

<i = P <i-1 + Q

Wd-1 = {fd e Vd : <fd, fd-1) = o V fd-1 e Vd-1} (6)

where (f, g) is the inner product of the functions f and g defined as

where P e R"iX"i-1 is a knot insertion matrix that can be computed through the Oslo algorithm [6], while Q e RK.'X(«.'-«.'-i) is a matrix that link the functions in the space V, to the B-wavelets function in W;_i, it is possible to find the procedure to compute the matrix in [1]. The matrix Q is computed such that the functions in are orthogonal to the functions in that is

($i-1,j,ii-1,k) = 0, V j = 1,...,ni-1, k = 1,...,ni - ni-1 < f, g) = f(t) g(t) dt, (7) (13)

It is possible to reconstruct the B-splines curve at level i as the sum of the B-splines and the B-wavelets curve belonging to the space i - 1

while it is the following weighted product for the decomposition of the NURBS function

<f, g)a = f a>2(t) f (t) g(t) dt, (8)

ri(u) = $(u) Ci = (<f>T-i(u) PT + fT-1 («) QT) Ci. (14) If the surface is described as the vector product of two B-

spline basis belonging to two spaces VUii and Vv>i, the decomposition of the tensor product space is [1]

Vu,i x Vv,i =(Vu,i-i x Vv,i-i) ® (V„,i-1 x Wvi-i)® (Vv,i-1 X Wu—) ® (Wuj-t x Wvj-t).

If we rearrange the control points of each dimension in a matrix of size n x m, control points matrices of the level i can be computed as

Ci,j =PuCi-i,j Pv + Pu DijQv + Qu D2,i Pu + QuD3,jQv, j = 1,...,r

where C,j is the matrix of control points if the decomposed surface is a NURBS. See [1] for more information on how to compute the matrices Ci_1>j, _D1>j, D2,j, D3,j.

3. Decomposition algorithm

The purpose of this paper is to decompose the original B-splines curve or surface at different levels of approximation and use these information to perform the sampling. In order to project a B-splines form the space V, into the sub-space Vi-1 we have to choose which interior knot, or knots, must be removed. In [7] a simplification algorithm has been proposed to incrementally decompose a B-spline curve. The initial step is to rank the knots according to their importance in the representation of the curve, this value is computed through the supremum norm

6j = |k(f) - rj-i(f)ll

where r^t) is the original curve and rj_ 1(i) is the curve without the j-th internal knots. The knots are then sorted according to the values 6j and the knot with the minimum value, or knots whose distances from the minimum are smaller than a threshold e, are removed. This process can be iterated until there are no more internal knots. Since are removed the less important values the details of the curve are described by the first B-wavelets curves.

A measurement method of the knot importance is also proposed in [8], the authors proposed to compute the importance of each knot as

6j = max {dist (Cj - {PC-1}j)}

where Сг-1 the matrix of the control points of the B-splines after removing the knot j, •j is the j-th row of the matrix • and dist(«) is a distance function [9]. This value compute the distance between the control points at different approximation levels, it describe the curve using only its control polygon. In this paper we use the value in Equation (18) with the euclidean norm as distance to compute the importance of each knots because it is faster compared to the computation of the supremum norm.

An example of a B-splines decomposition is shown in Figure 1. The blue curve is the original curve and the orange curves represent the simplified curves; after each iteration of the B-wavelets decomposition some details of the curve are lost. The last B-splines curve is a three degree curve without internal knots.

(a) (b)

Fig. 1: (a) first decomposition level; (b) last decomposition level

The maximum number of the decomposition levels (functions) can be checked through the number of the internal knots of the NURBS parameter space. It is possible to increase the number of B-wavelets basis functions modifying the number of the knots with a refinement operation [5].

A NURBS surface is decomposed ranking pair of knots in u and v direction. For each couple of knots the rank is based on [8]

6ki = max {dist Ci - {Pu,Ci-1,jPTv,i(19)

where Ci-1, j the matrix of the control points of the NURBS after removing the j-th knot, *k'1 is the matrix entry in position kl and dist(«) is a distance function, as for in the curve decomposition we use the euclidean distance. An example of a surface [10] decomposition is shown in Figure 2, where there is the initial surface and the B-splines surface after the first iteration, the colour represents the value of the B-wavelets surface.

Fig. 2: Surface decomposition example

4. Sampling strategy

In this Section we will explain the sampling strategy proposed in this paper. In order to take into account both the length of the curve and its complexity, we propose a sampling method based on a curve parameterisation computed as a combination between the arc length and the B-wavelets basis.

The arc length parameterisation of a parametric curve r(t),

with t e [i1, tn+1] can be computed as [11]

l(t) = ||r'(u)|| du.

Since the difference between successive curve approximation is measured by the value of the B-wavelets curve, we use the square of the B-wavelets basis functions as measurement of the complexity of the shape. The sum of the values at the i-th level

possible to observe that there are more points where there are abrupt changes of the curve and on the two boundaries, while on the flat zones the spacing between the samples increases. The proposed surface sampling, as for the curve one, is based on a re-parameterisation of the surface based on a combination of its area and square of the value of the B-wavelets basis. The area of a parametric surface can be computed as [11]

n Un+1 nv,

Ju1 Jv1

||ru(u, v) X rv(u, v)|| dudv

E.(t) = £ <¿2(t).

where x is the cross product.

In order to measure the complexity of the surface, we compute the sum of the square value of the B-wavelets basis functions at the i-th level

Since this value is a function on the B-splines curve r(t), the cumulative integral of E.(t) is computed

Ec(t) = f E.(t) Hr'(u)H du. (22)

n.+1 _n. m.+1 _m.

E.(u, v) = ^ ^ <2(u, v)<2(u, v),

j=1 k=1

so the cumulative integrals of the i-th B-wavelets basis functions along the surface r(u, v) are

This value represents the energy of of the B-wavelets basis function along the curve. The values of the cumulative integrals must be merged in a function to represent the complexity of the decomposed curve; we compute this "mean" function as

E(t) = - £

n ¿—t

¡=1 E\(tn+1)

The value of the i-th B-wavelets is divided by its maximum to assign the same importance to each approximation scale. Once both the cumulative arc length and the B-wavelets based function have been computed, the parameterisation based on a mixture between these two quantities can be computed as

ECu(u) =

fu rvm+\

E.(s, v) ||rs(s, v) X rv(s, v)|| dv (27)

i,u^"> ~ I I '/ll'sVJ> ")X 'v\

Ju1 Jv1

Jrv fun+1

ds E.(u, s) ||ru(u, s) X rs(u, s)|| du. (28)

We can compute the marginal cumulative areas, Au(u) and Av(v), as the previous equations without the weights E.(s, v). The function representing the marginal complexity of the surface are computed as

&(.) = 1 £

El.(') 1 max^ Ec •(•)

-=2 («£>+e«

where there is a normalisation to sum quantities of different approximation levels. We finally compute the marginal parame-terisations as

where the arc length is coded in the range [0, 1] because we are adding quantities with different measurement units, and also to assign equal weight to both the values.

An example of the uniform sampling of 60 points with the parameterisation in Equation (24) is shown in Figure 3. It is

1 / Au(u) Pu(u) -- + Eu(u)

2 \Au(um+1 )

PW> = aAWD + Ev(v».

Fig. 3: B-wavelets' sampling (60 points)

where the areas are divided by their maximum values to give to both quantities the same importance. Once we have computed these two functions, a sampling method can be computed through the vector product of two uniform sampling along the parameterisations pu(u) and pv(v). An example of the proposed sampling is shown in Figure 4. This sampling is uniformly distributed in the internal part of the surface, while the density of the points is higher close to the boundary.

5. Reconstruction results

In this Section we evaluate the performances of the sampling strategy described in Section 4. The reconstruction of the profile in Figure 3 and a profile based on the simulation of a manufactured surface are firstly computed.

Fig. 4: B-wavelets surface sampling (20 x 20 points)

Then the performances of the proposed sampling method is evaluated with the reconstruction of three free form surfaces.

We use the Root Mean Square Error (RMSE) as performance indicator

££M-

i=i j=i

j(s))2

Fig. 5: RMSEs as a function of the number of the samples

Fig. 6: (a) B-wavelets sampling of a simulated milled profile (100 points), all the values are in mm; (b) RMSEs of the compared sampling strategies

Figure 6b. The values of the arc length and the uniform param-eterisation are very close, while the reconstruction with LHS sampling sometimes achieve better results compared to the previous two. The proposed sampling outperforms other methods with all the analysed sample sizes because it is able to add more points where there are rapidly changes of the curve.

where n is the number of the points, d is the dimension of the B-spline curve or surface, rjs) is the true, known, value and Tij(s) is the estimated point with the reconstruction model.

5.1. Curve reconstruction

In the curve reconstruction test cases, we analyse four different sampling strategies:

• uniform sampling on arc length parameterisation, Arc length;

• uniform sampling according to the parameterisation in Equation (24), B-wavelet;

• uniform on the parameter domain t, Uniform;

• random latin hypercue sampling [12], LHS.

The reconstruction algorithm used in this paper is the piece-wise cubic spline interpolation implemented in MATLAB® [13] through the cubucinterp function. In this analysis we assume that the true value of the parameter of the curve is known, i.e. we don't consider the point cloud parameterisation.

The values of the RMSEs as a function of the number of samples of the curve in Figure 3 are shown in Figure 5; all the methods, except the one based on the LHS sampling, achieve comparable results. We now present a comparison based on the reconstruction of a simulated manufactured profiles [14]. Starting form the measured points, we have reconstructed the curve with a free knots spline regression, where the position of the knots is set according to the changes of the profile [15].

The reconstructed profile with 100 sampled points is plotted in Figure 6a, and the corresponding RMSEs values are shown in

5.2. Surface reconstruction

The proposed sampling method is tested with the reconstruction of three surfaces, the free form surface in Figure 4, a structured surface and a surface with grooves where there are abrupt changes along one coordinate. We analyse six different sampling strategies:

• uniform sampling based on the cumulative marginal areas, Area;

• uniform sampling based on the marginal parameterisations in Equations (30-31), B-wavelets;

• uniform sampling on the parameters space u-v, Unif;

• Latin Hypercube Sampling, LHS;

• Hammersley sampling [16], Hammersley;

• Halton sampling [16], Halton.

Since we assume deterministic surfaces, i.e. without measurement error, the reconstruction is performed through the piece-wise cubic interpolation implemented in the fit function of MATLAB1®.

The RMSEs of the reconstruction of the first surface are shown in Figure 7; the area based sampling achieve better results when the sample size is small, while the proposed sampling method has similar results when the number of samples is equal or greater than 225 points. The proposed method has not the same results of the area based method because the parame-terisation adds more points close to the boundary of the domain and when the sample size is small the reconstruction is good in those parts, but the error is bigger elsewhere.

In Figure 8a is shown the second surface analysed along with

a sampling of 40 x 40 points and its reconstruction. It is a C2 structured surface with a small radius of curvature on the edges.

ID Д X 40 И

Fig. 7: RMSEs of the compared sampling strategies

The corresponding RMSEs are shown in Figure 8b. In this second test case the sampling based on the B-wavelet gets the best performance, while the area based sampling has worst results also compared to the Unif and Hammersley samplings.

(a) (b)

6. Conclusion

In this paper we have proposed a sampling method based on a parameterisation that is the mean between the curvilinear abscissa, or area in the surface case, and the square of the values of the B-wavelets basis functions. This parameterisation has been tested with other common sampling methods and in different test cases. It has been proved that it is able to achieve better prediction results if the reconstructed shape is a profile. If the shape to be reconstructed is a surface, the proposed sampling method achieves better results if the analysed surface has abrupt changes. If the analysed surface has slow changes, the area based parameterisation ensures a better reconstruction error if the sample size is small. In this paper we have taken into account the value of the basis function, but it is also possible to build a parameterisation based on the differential property, such as the curvature of the B-wavelets curves or surfaces at different levels.

Acknowledgements

The authors gratefully acknowledge the UKs Engineering and Physical Sciences Research Council (EPSRC) founding of the EPSRC Fellowship in Manufacturing: Controlling Variability of Products for Manufacturing (Ref:EP/K037374/1).

Fig. 8: (a) B-wavelets sampling of a structured surface (40 x 40 points); (b) RMSEs of the compared sampling strategies

The third surface analysed, see Figure 9a, is a surface with three grooves, it has some abrupt changes of the surface along the x direction. Since the complexity of the surface, that is described by the internal knots, is different in the two directions, we propose to use

References

■ (ntot - 4)

where nu and nv are the number of points in u and v directions, IXI is the closest integer of x, nk is the number of the internal knots in direction i and ntot = nu + nv is the total sample size in the marginal directions. The RMSEs of this test case are shown in Figure 9b, there are only the values of the Area, B-wavelts and Unif sampling because the others are too high. Also in this reconstruction scenario the proposed method achieves better prediction results.

Fig. 9: (a) B-wavelets sampling of surface with grooves (100 X 2 points); (b) RMSEs of the compared sampling strategies

[13 [14

Lyche, T., Mrken, K., Quak, E.. Theory and algorithms for non-uniform spline wavelet. In: Dyn, N., Leviatan, D., Levin, D., Pinkus, A., editors. Multivariate Approximation and Application. Cambridge University Press; 2001, p. 152187.

Chui, C., Lian, J.. Multilevel structure of nurbs and formulation of nur-blets. In: Wavelet Analysis: Twenty Years' Developments. 2002,. Cottrell, J., Hughes, T., Bazilevs, Y.. Isogeometric Analysis: Toward Integration of CAD and FEA. 2009.

Petkovski, M., Bogdanova, S., Bogdanov, M.. A simple adaptive sampling algorithm. In: 14th Telecomunications forum TELFOR. 1006, p. 329-332.

Piegl, L., Tiller, W.. The NURBS Book. Springer-Verlag New York, Inc.; 1997.

Lyche, T., M0rken, K.. Making the oslo algorithm more efficient. SIAM Journal on Numerical Analysis 1986;23(3):663-675. Lyche, T., M0rken, K.. A data reduction strategy for splines. IMA Journal of Numerical Analysis 1988;8:185-208.

Wang, W., Zhang, Y.. Wavelets-based nurbs simplification and fairing. Computer Methods in Applied Mechanics and Engineering 2010;199(5-

8):290300.

Arkhangel'skii, A., Pontryagin, L.. General Topology I: Basic Concepts

and Constructions Dimension Theory. Springer; 1990.

X3D, . Example archives. 2015. URL: http://www.web3d.org/x3d/

content/examples/Basic/NURBS/.

do Carmo, M.. Differential Geometry of Curves and Surfaces. Prentice-Hall; 1976.

McKay, M., Deckman, R., Conover, W.. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 1979;21(2):239-245. MATLAB. Natick, Massachusetts: The MathWorks Inc.; 2014. NPL, . Softgauges. 2015. URL: http://resource.npl.co.uk/ softgauges/AboutSG.htm.

Luong, B.. Free-knot spline approximation. 2009. URL:

http://it.mathworks.com/matlabcentral/fileexchange/

25872-free-knot-spline-approximation.

Wong, T.T., Luk, W.S., Heng, P.A.. Sampling with hammersley and halton points. J Graph Tools 1997;2(2):9-24.

n;Ê + n