Improving the approximation quality of tensor product B-spline surfaces by local parameterization

2.1 B-spline curves

B-spline curves are parametric curves of the form [4]:

According to Equation (1), each coordinate component is an explicit function of the curve parameter u, which locates a point C(u) on the curve. The one-dimensional parameter space is mapped into the Cartesian space by [4]:

Thus, a curve point C(u) is computed as the weighted average of n _P + 1 control points P _i , i = 0, …, n _P . The corresponding weights are defined by the B-spline basis functions N_i,p(u), which can be computed recursively by means of the Cox-de Boor-algorithm (see [25, 26]). The B-spline basis functions are completely defined by their degree p and a non-decreasing sequence of r + 1 knots u k * , typically aggregated in a non-periodic knot vector U [4]:

This knot vector subdivides the parameter space into knot spans. Due to this subdivision, the basis functions as well as the B-spline curve itself are piece-wise polynomial curves, which are assembled at the knots.

2.2 Tensor product B-spline surfaces

B-spline surfaces are parametric surfaces and, thus, coordinate-wise defined functions of the two surface parameters u and v [4]:

In analogy to Equation (2), the parameter space is mapped into the three-dimensional Cartesian space by [4]:

Similar to the curves, a surface point S(u, v) is computed as the weighted average of a bidirectional net of (n _P + 1) ⋅ (m _P + 1) control points P_i,j. The respective weights are defined by bivariate basis functions, which result from the multiplication of two one-dimensional basis functions N_i,p(u) and N_j,q(v) with degrees p and q, respectively. The basis functions are defined on the knot vectors U and V, respectively [4].

The B-spline surface’s representation in Equation (5) is a TP representation. This form of B-spline surfaces represents them as an infinite number of B-splines curves that run into two different directions [4]. Curves that lie on this surface and run either in u- or in v-direction are called isoparametric curves. They are constructed by fixing either u or v to be a constant. Due to the strong relationship between B-spline curves and surfaces, the properties of B-spline curves can be transferred to a large extent to TP B-spline surfaces.

2.3 Skinned B-spline surfaces

The above described relationship between B-spline curves and TP B-spline surfaces is especially useful when constructing certain kinds of B-spline surfaces, one of them being the skinned B-spline surface. In this context, “skinning” or “blending” is the process of assembling B-spline curves so that they form a B-spline surface [4]. The curves used for the skinning process are then isoparametric curves on the resulting surface.

Starting point for the construction of skinned B-spline surfaces is a set of m _P + 1 section curves C _j (u) (j = 0, …, m _P ), each of the curves being available in the form of Equation (2) [4]. Usually, these curves are planar cross sections (see Figure 1). In accordance with their denotation, they are defined in u-direction of the surface to be constructed. As a prerequisite for the skinning process, all section curves have to be defined on the same knot vector U and are required to have the same degree p. If these conditions are not met, the curves can be adjusted by means of knot refinement and degree elevation (cf. e.g. [4]). Both algorithms allow for the adaption of the respective parameter group without changing a curve’s shape (see Figure 2).

Figure 1:

Section curves for a skinned surface to be constructed (coloured curves: B-spline curves, black crosses: control points).

Figure 2:

Section curves for a skinned surface to be constructed after knot refinement and degree elevation (coloured curves: B-spline curves, black crosses: control points).

With these section curves being isoparametric curves on the TP B-spline surface S(u, v) to be constructed, they can be related to the B-spline surface by fixing the surface parameter v to be a constant v _k [27]:

Comparing Equation (7) with Equation (2) reveals strong similarities between the two equations: The control points in Equation (2) are control curves C _i (v _k ) in Equation (7):

Fixing the value v = v _k , thus, results in n _P + 1 points on these control curves. These n _P + 1 points are the control points of the isoparametric section curve S(u, v _k ) [27].

In order to construct control curves based on given section curves, three parameter groups have to be specified: The degree q, the knot vector V and the control points P_i,0, …, P i , m P . The degree has to meet the requirement q ≤ m _P and can be chosen arbitrarily apart from that [4]. Equally spaced knots are an appropriate choice when using equally spaced section curves [4], which leaves the determination of the control points. As is apparent from Equations (7) and (8), the control curves interpolate the control points of the section curves. Thus, the control points can be determined by solving this interpolation problem [27]. Figure 3 presents the control curves (grey curves) constructed in this manner for the example section curves.

Figure 3:

Control curves (grey curves) constructed by means of the coloured section curves.

Following Equation (7), the control mesh of the skinned surface of degrees p and q, which is defined on the knot vectors U and V, is then obtained by assembling the control curves’ control polygons to a (n _P + 1) ⋅ (m _P + 1) control grid [27] (see Figure 4 for the final surface).

Figure 4:

Skinned surface and the section curves used for construction.

2.4 Point cloud modelling with B-spline curves and surfaces

In the context of geodetic laser scanning, the fitting of B-spline curves and surfaces through a set of three-dimensional data points X k = [ x k , y k , z k ] T (k = 1, …, n _l ) is of major interest. Usually, only the control points’ three-dimensional positions are estimated [4, 5]. This strategy results in a linear relationship between the observations l = […, X _k , …] and the unknown control points, as exemplified for a B-spline curve in Equation (9). (The functional model of the curve approximation can be directly transferred to the estimation of B-spline surfaces.)

In case all remaining parameter groups (degree(s) p (and q), number of control points (n _P + 1) (⋅(m _P + 1)), knot vector(s) U (and V) as well as the curve/surface parameters u _k (and v _k )) are known, the control points can then be estimated in a linear Gauß–Markov model (cf. e.g. [28]).

However, usually, there is no a-priori knowledge regarding the remaining parameter groups, making a multi-step procedure necessary [5]: Initially, appropriate curve/surface parameters are allocated to the observed data points. These curve/surface parameters provide the basis for the computation of the knot vector(s) (see [4] for possibilities to determine the knot vector(s)). The determination of the optimal number of control points is a model selection problem and can be solved e.g. by means of information criteria or by means of structural risk minimization [29, 30]. The latter also allows the degree to be taken into account in model selection [31]. Otherwise, the use of cubic B-splines with p = 3 (and q = 3) is a generally accepted choice [4].

Especially the curve/surface parameterization has turned out to be a very critical point within B-spline-based point cloud modelling, as the curve/surface parameters locate the observations on the curve/surface to be estimated. Thus, the parameterization requires prior knowledge about the curve/surface to be estimated. The better the prior knowledge is, the better is also the point cloud approximation.

In spite of these difficulties, a variety of parameterization strategies exists. Uniform, chordal and centripetal parameterization can be seen as standard strategies for ordered point clouds [32]. As laser scanning profiles can usually be ordered, these strategies can be used when modelling point cloud profiles by means of B-spline curves. However, for areal laser scanning point clouds, an ordering of the data points can usually not be realized. Hence, base surface parameterization is the standard approach for unordered areal point clouds [32]. Starting with a parametric base surface that roughly approximates the point cloud to be parameterized, surface parameters are obtained by projecting the point cloud onto the base surface. Coons Patch, which is constructed of the point cloud’s boundary curves, has been proven to be an appropriate initial base surface for parameterizing point clouds representing a variety of engineering geodetic objects [33]. In case Coons Patch does not satisfactorily approximate the point cloud, an iterative procedure, which uses the resulting estimated and reparameterized B-spline surface as a new base surface, can considerably improve the parameterization. Usually, however, boundary constraints have to be introduced in order to prevent the surface’s degeneration during the iteration [33].

4.1 Test specimen

A test specimen representing a B-spline surface with 9 ⋅ 7 control points is available at TU Wien [29] (see Figure 6). Due to the highly accurate known shape of the specimen, it is perfectly suitable for evaluating newly developed modelling approaches. Simulated point clouds representing this test specimen are used in this paper to evaluate the developed approach. The advantage of simulated data over measured data is manyfold: No systematic effects e.g. caused by the strongly varying incidence angles falsify the evaluation, there is no influence of an incomplete stochastic model on the estimation results and a ground truth, against which the results can be compared, is available.

Figure 6:

B-spline shaped test specimen [29].

As a perfect approximation of the test specimen can only be expected when the nominal surface parameters are used, an unnoisy point cloud of the test specimen is used as a first data set (see Figure 7). In this way, a comparison that is undistorted by noise can be made with the nominal surface. This data set Specimen _u is complemented by a noisy data set Specimen _n , which is generated by adding white noise (σ _x = σ _y = σ _z = 1 mm) on the simulated point cloud. With this publication’s focus being on the establishment of an improved functional model rather than on the investigation of the stochastic model, this is an appropriate simplification.

Figure 7:

Non-noisy point cloud of the test specimen (blue) and rectangular segment (red) used for the investigations within this publication.

In addition to the surface’s undulation, the specimen is a challenging object due to the non-rectangular boundary. In the investigations conducted in this paper, the handling of these types of boundaries is not focused on. For this reason, not the entire point cloud of the specimen, but a rectangular segment (red point cloud in Figure 7) is used for the initial developments. The approach will then be extended to non-rectangular boundaries in future.

Both point clouds represent an undulating freeform surface. The point clouds themselves are free of outliers and data gaps and are characterized by a very regular distribution of the data points.

4.2 Beach scene [34]

The third point cloud is part of a published laser scan data set of the Kijkduin beach-dune system in The Netherlands [34]. The data set originates from a long-term monitoring of a sandy beach at the North Sea coast in The Hague, The Netherlands, by means of permanent laser scanning. For this purpose, a long-range laser scanner has been installed on top of a hotel at about 38 m height above mean sea level so that it overlooks the beach and the dunes. During a 6-months period in 2016/2017, hourly scans have been conducted. The resulting point clouds are georeferenced and checked against RTK-GNSS and airborne laser scanning (see [34] for further details regarding the data acquisition and preprocessing). From the resulting 4082 point clouds, one has been arbitrarily chosen (acquired on 11.11.2016) for the investigations conducted within this paper. In order to obtain a point cloud with clearly defined boundaries, a rectangular part was cut out of this point cloud (see Figure 8). This step simplifies the initial investigations carried out in this publication and will be omitted in future investigations.

Figure 8:

Preprocessed point cloud of the beach scene (data source: [34]).

The resulting point cloud is characterized by a lot of local structures (the dunes) and in some regions by irregular point densities and even data gaps. Both properties are caused by the scanner’s flat viewing angle and the resulting occlusions.

4.3 Soil scene

The fourth point cloud originates from a long-term soil erosion monitoring programme, which aims to detect small scale erosion under field conditions in northern Germany. For this purpose, small plots with dimensions of 2 × 3 m are installed on managed cropland (maize) and marked by means of ground sleeves [35]. Starting with the scans immediately after seeding, the plots have been acquired on a weekly basis (11th May to 8th June 2022), resulting in point clouds of five measuring epochs. In each epoch, the plots are acquired by means of a TLS from two different viewpoints. In order to align the acquired point clouds into a joint coordinate system, a target-based registration is performed. Finally, as the maize plants that have grown over time are not of interest for the actual analysis of the soil, a ground filtering is conducted (see [35] for more information regarding the data acquisition and preprocessing). One of the acquired point clouds serves as the data set Soil in this publication. This point cloud is presented in Figure 9.

Figure 9:

Preprocessed point cloud of the soil scene (data source: [35]).

The acquired point cloud is characterized by a lot of local structures (sewing rills and soil roughness). Furthermore, due to the removed vegetation and due to occlusions, small data gaps exist. Finally, as the ground filtering is not yet fully developed, a multitude of outliers that are caused by data points representing vegetation exist.

5.1 Basic ideas

TP B-spline surfaces as described in Section 2.2 can be interpreted as B-spline surfaces that are composed of an infinite number of B-spline curves that run into two different directions [4]. Each of these curves is an isoparametric curve. Usually, TP B-spline surfaces are attributed the weakness of having a global character which arises by the global definition of most of the parameter groups (e.g. the knot vectors, the degrees and the number of control points): Each of the isoparametric curves running in the same direction is defined on the same knot vector and has the same degree as well as the same number of control points. When determining the best-fitting TP B-spline surface as described in Section 2.4, these parameter groups are simultaneously determined for the entirety of isoparametric curves. This strategy considerably restricts the flexibility of TP B-spline surfaces.

Taking these considerations into account, the common strategy to determine best-fitting TP B-spline surfaces is adapted in the following paragraphs: Starting point is the determination of best-fitting section curves (Section 5.2). These section curves are afterwards used to construct skinned base surfaces (Section 5.3). Due to the independent determination of the isoparametric curves, the resulting skinned TP B-spline surfaces receive a local character and, thus, an increased flexibility. By using these skinned surfaces as base surfaces (Section 5.4), considerably improved surface parameters can be allocated to the observations, which implies a considerably improved point cloud approximation. Iterative reparameterization (cf. [33]), can further improve the point cloud approximation. The price paid for the base surfaces’ increased flexibility is a globally increased number of control points. However, as the skinned surfaces are only an intermediate result, this does not have an impact on the final TP B-spline representation.

Thus, the basic idea of this strategy is to store important information about the surface’s shape – especially about local variations – in the surface parameters. As a result, a comparatively simple and compact representation in form of a TP B-spline surface is possible. As base surfaces with a local character are used for this purpose, the developed B-spline surfaces are designated as locally parameterized (LP) B-spline surfaces.

The individual steps outlined above are schematically presented in Figures 10 and 11 and described in more detail in the following sections. It is worth noting that in all steps described below a simplified stochastic model is used that does not consider correlated observations, although some of the data sets do contain correlated noise (cf. Table 1). This simplification is made as a complete stochastic model is not available for the measured data sets, and as the focus of this publication is on the functional rather than on the stochastic model. A complete stochastic model will be incorporated in future research.

Figure 10:

Flowchart for the estimation of locally parameterized B-spline surfaces. The red coloured parts are presented in detail in the flowchart in Figure 11.

Figure 11:

Flowchart for the estimation of section curves. (PC: point cloud).

5.2 Estimation of section curves

In order to determine the section curves, the three-dimensional point cloud is initially subdivided into n _s point cloud stripes. Since all data sets used in these investigations are oriented parallel to the x/y plane, this subdivision can be realized via the x- or y-coordinates of the data points. Optimally, the direction is chosen in a way that the stripes run orthogonally to the existing surface relief. Further impacts of the initial subdivision will be regarded in Section 5.4.

The points of each stripe are then projected into a mean vertical plane, resulting in n _s two-dimensional point cloud profiles. A reasonable number of profiles n _s has to be chosen. The larger n _s , the narrower the point cloud stripes and the fewer points are contained in the resulting profiles. If the point cloud stripes are chosen too narrow, there are not enough points in each profile. Contrarily, if n _s is chosen too small, the profile averages over too wide a range of the point cloud. In both cases, the point cloud cannot be appropriately approximated by the section curves. In the initial investigations conducted within this publication, the number of profiles n _s is chosen manually. In future, however, this step is going to be automized by taking the considerations above into account.

Afterwards, a best-fitting B-spline curve is determined for each of the profiles. Within the scope of the investigations, different approaches for the curve estimation, which are introduced in the following paragraphs, are combined.

5.2.1 Centripetal curve parameterization

As a starting point, centripetal parameterization [32] is used for assigning curve parameters to the data points. Based on these curve parameters, an initial B-spline curve is estimated as described in Section 2.4.

In Figure 12, one of the point cloud profiles of the data set Beach, which is characterized by a lot of undulations, can be seen as an example (blue dots). Two centripetally parameterized B-spline curves are presented as examples in addition to the point cloud profile: The red curve with 12 control points is too simple to appropriately represent the point cloud’s undulations. An increase of the number of control points, which increases the curve’s flexibility, clearly improves the approximation quality (yellow curve with 22 control points). From a pure visual point of view, centripetal parameterization, thus, yields satisfying results. However, a closer look at the adjustment’s residuals reveals an unexpectedly large sum of squared residuals ( ɛ ^T ɛ = 296.521 m²). Comparing this value with the sum of squared Euclidean distances between data points and curve (d ^T d = 5.342 m²), a large discrepancy becomes visible. This discrepancy can be explained by the way of computing the residuals: The point correspondences between observations and estimated curve points, which are needed to compute the residuals, are defined by the curve parameter u. The residuals, therefore, are not necessarily perpendicular to the curve. By contrast, the curve parameter is not taken into account, when computing the Euclidean distances; rather, for each data point the nearest neighbour on the curve is looked for. However, with the sum of squared residuals being the target value to be minimized during optimization, it should be striven for residuals that are perpendicular to the curve and, hence, correspond to the Euclidean distances. The large discrepancy between ɛ ^T ɛ and d ^T d, therefore, suggests a poor parameterization.

Figure 12:

Results of the curve estimation when using centripetal parameterization. The point cloud profile (blue dots, data set Beach) is approximated by means of B-spline curves with 12 (red curve) and 22 (yellow curve) control points.

5.2.2 Iterative curve reparameterization

Adapting the idea of Ref. [33] to curve parameterization, an iterative reparameterization approach is developed to improve the curve parameterization. This step is necessary, as the initial curve parameters are obtained by means of a centripetal parameterization using noisy observations. As a consequence, the corresponding parameter space does not appropriately reflect the actual shape of the curve. By iteratively allocating new curve parameters, the parameter space is gradually adapted to the actual shape of the curve.

The procedure starts with a B-spline curve received by means of centripetal parameterization. Afterwards, curve points lying on the estimated B-spline curve are generated with a high sampling rate. By centripetally reparameterizing the resulting curve points, the parameter space of the estimated curve is adjusted to the curve’s actual shape. Afterwards, the point cloud is projected onto this curve, and improved curve parameters that represent the new parameter space are allocated to the observations. This step is realized by finding the two closest curve points for each data point and by linearly interpolating the curve parameter to the data point. The new parameter values are used in the subsequent iteration step to estimate improved control points and, thus, an improved B-spline curve. This curve is then reparameterized as described above and is then used as a new base curve. This process is carried out until either a predefined number of iterations is reached or the parameter values do not change any longer.

From a pure visual point of view, the shapes of the resulting B-spline curves do not considerably change compared to the classical centripetally parameterized B-spline curves. The difference between the results becomes only apparent when comparing the curve parameters, as is done in Figure 13. For the example profile presented in Figure 12, the maximal absolute difference in the resulting curve parameters d _u = |u_centri − u_iter| is up to 0.025.

Figure 13:

Influence of the parameterization strategy on the resulting curve parameters. The estimated curve points are coloured according to the difference d _u = u_centri − u_iter.

These differences in the allocated curve parameters have a strong influence on the sum of squared residuals, as is apparent from Table 2, listing ɛ ^T ɛ and d ^T d for ten example profiles of the data set Beach using both, centripetal parameterization and iterative reparameterization. For all ten profiles, the sum of squared residuals can be clearly reduced by iteratively adapting the parameter space to the actual shape of the curve. As a consequence, the sum of squared residuals approaches the sum of squared distances, indicating residuals which are considerably closer to being perpendicular to the curve than the residuals resulting from centripetal parameterization. The improved parameterization also improves the approximation of the point cloud in terms of a smaller sum of squared distances.

Table 2:

Sum of squared residuals ɛ ^T ɛ and sum of squared Euclidean distances d ^T d for ten example profiles of the data set Beach (centripetal parameterization: _centri and iterative reparameterization: _iter).

Profile	( ε T ε ) centri [ m 2 ]	( d T d ) centri [ m 2 ]	( ε T ε ) iter [ m 2 ]	( d T d ) iter [ m 2 ]
1	157.4904	1.8668	2.5762	1.6805
2	144.1247	1.7005	2.4298	1.5273
3	148.0414	1.4369	2.3922	1.3891
4	136.1523	1.5947	2.5074	1.4449
5	141.6964	2.3913	2.9512	1.8398
6	138.1903	2.2434	3.1654	2.1629
7	151.6335	2.9497	3.7015	2.5682
8	165.7931	3.5175	4.3483	3.3354
9	158.6155	3.8832	4.7647	3.8557
10	175.6595	5.1642	6.0383	5.0038

5.2.3 Robust curve estimation

Laser scan data representing natural surfaces often contain a considerable amount of outliers. In case of the data set Soil, for example, not only the soil surface is acquired by the laser scanner, but also the plants growing on it. As only the soil surface is of interest in the corresponding investigations, a ground filtering is conducted (see [35] for further information). However, even after this step there are usually still unwanted vegetation points in the point cloud, which are outliers with respect to the subsequent investigations. To make point cloud modelling with B-spline curves and, thus, the parameterization, more robust with respect to outliers, a robust procedure is therefore used in the presented investigations. For this purpose, the approach proposed by Ref. [36] and briefly described below is used in the presented studies.

As described in Section 2.4, the control points of a B-spline curve are estimated in a linear Gauß Markov model when fitting a B-spline curve to a point cloud profile. Within the estimation procedure, the data points’ uncertainty is taken into account by means of the observations’ variance covariance matrix (VCM) Σ ll = σ 0 2 ⋅ Q ll , with Q_ll being the corresponding cofactor matrix and σ 0 2 being the a-priori variance factor. In order to detect outliers in the data, Ref. [36] proposes a step-wise approach: The approach’s basic idea is that observations with larger residuals ɛ _i should have a smaller influence on the estimated curve. The influence function of the Huber estimator used within these investigations is given by [36]:

(10) ψ ( ε i ) = ε i sign ( ε i ) ⋅ c for | ε i | < c for | ε i | ≥ c .

According to Equation (10), observations with residuals that are smaller than a tuning constant c have the same influence as they have within classical least squares estimation. In case a residual’s magnitude exceeds this constant, the corresponding influence is limited by c [36]. The robust estimation of the B-spline curve is then implemented as described in Ref. [37].

Starting point of the procedure is an initial least squares estimation without consideration of any weights. The control points of the B-spline curve, summarized within the vector of unknowns ϑ , the observations’ residuals ɛ as well as the variance factor σ as median absolute deviation (mad) are estimated on the basis of the design matrix A and the observations l [36]:

ϑ ̂ ( 0 ) = A T A − 1 A T l ε ( 0 ) = A ϑ ̂ ( 0 ) − l σ ( 0 ) = 1,4785 ⋅ mad ( ε ) .

In the next step, a weighted estimation is used to iteratively determine improved values for ϑ ̂ , ɛ and σ (cf. Algorithm 1). The outer while loop of this iteration is interrupted as soon as the residuals do not longer change considerably from one iteration step to another. The iteratively computed weights p _i used to establish the weight matrix Q ll − 1 take into account the influence function ψ as a function of the residuals and the variance factor (cf. Algorithm 1).

For data points with for |ɛ _i |≥ c, the influence function fulfils

ψ v i ( k ) σ ( k ) < v i ( k ) σ ( k ) ,

resulting in weights p _i < 1. The remaining data points are fully taken into account with p _i = 1.

In principle, this strategy is very similar to the commonly used outlier detection based on standardized residuals within geodetic networks, where the tuning constant c is chosen according to the observations’ probability distribution.

Algorithm 1

While-loop of iterative reweighted least squares algorithm [36].

The influence of the robust curve estimation is presented as an example in Figure 14. The reduced influence of data points with large residuals (indicated by the black rectangles in Figure 14, right) is clearly visible in a changed shape of the estimated B-spline curve compared to the result obtained by classical least squares estimation (Figure 14, left). Despite of the successful application of the robust approach, Figure 14 simultaneously highlights two challenges of the robust curve estimation: At first, there is a large influence of the tuning constant c that directly determines the number of detected outliers. In the example, the tuning constant has been chosen very pessimistically, so that the majority of data points is downweighted. A careful choice of c, taking the stochastic properties of the scanner used as well as possible systematic effects into account, is therefore necessary. Secondly, Figure 14 (right) highlights the difficulty to distinguish between local variations of the surface and outliers, which becomes especially challenging when not only single outliers but accumulations of outliers occur in the point cloud. This differentiation will be part of future research.

Figure 14:

Comparison between classical and robust B-spline curve estimation. Data points with large residuals (marked by the squares) are downweighted in the robust approach (right).

Figure 11 summarizes the entire strategy for the estimation of section curves.

5.3 Construction of skinned base surfaces

The section curves determined as described in Section 5.2 are then used to construct skinned surfaces (cf. Section 2.3). Starting with centripetally parameterized section curves (cf. subsection 5.2.1), the resulting skinned surfaces show artefacts in some regions (see Figure 15) and are not at all usable for further processing steps. These artefacts emerge, as the control points of the control curves are linked to the parameter space of the respective section curve via the corresponding knot spans. If the section curves’ parameter spaces considerably differ from each other, the control points of the control curves jump from one section curve to another in direction of the skinning parameter. This behaviour results in a surface with artefacts.

Figure 15:

Skinned surface for the data set Beach based on centripetally parameterized section curves and emerging artefacts (encircled in red).

The artefacts can already be considerably reduced by using iteratively reparameterized section curves (cf. Section 5.2.2). However, also in this case, the estimation of the section curves is independently realized. On the one hand, this is exactly the motivation for using skinned surfaces, as this strategy gives skinned surfaces their local character. On the other hand, this neglects the fact that neighbouring curves cannot be completely independent of each other, provided that continuous surfaces are assumed. Hence, the independent curve parameterization yields discontinuities in the resulting skinned surface’s parameter space. Taking the considerations above into account, an approach is developed in which the neighbouring section curve is involved when estimating a section curve.

The procedure starts with the estimation of one of the section curves, which can be arbitrarily chosen, as initial reference curve C ref 0 ( u ) . In the presented case, the middle curve is used, but alternatives are conceivable. When estimating the reference curve, it is important that it correctly describes the point cloud profile and, for example, that outliers are ignored or downweighted. For this reason, a robust estimation (cf. Section 5.2.3) should be conducted.

Starting with the profile nearest to the reference profile, afterwards, the remaining n _s − 1 point cloud profiles are passed through one after the other until the point cloud’s border is reached. Each point cloud profile that is regarded during this iteration is then approximated by means of a B-spline curve. During this step, the corresponding and already estimated neighbouring curve serves as base curve for parameterization (cf. Figure 11) and is involved in the respective curve estimation as described in the following paragraph.

For estimating a section curve C _i (u) by taking into account the neighbouring curve C_i−1(u), the latter is shifted to the ith profile in the first step. The shifted neighbouring curve is then sampled with a high resolution. Analogously to the iterative reparameterization procedure described in Section 5.2.2, initial curve parameters for the ith point cloud profile are then generated by projecting each point of this profile onto the shifted neighbouring curve. Hence, rather than using the centripetal parameterization as a starting point for the iterative reparameterization, the parameterization of the neighbouring curve is used, which creates a dependency between these two curves. By passing through the profiles step by step, this dependency propagates from one curve to the next.

After having allocated initial curve parameters to the ith point cloud profile, the curve parameters can be further improved by iteratively reparameterizing according to Section 5.2.2. The iteration can be interrupted after n_iter steps in order to maintain the dependency between the two adjacent curves. The smaller n_iter is chosen, the stronger is the influence of the neighbouring curve. Contrarily, when using a large n_iter, the influence of the adjacent curve decreases.

Figure 16 presents the resulting skinned surface for the data set Beach when interrupting the iteration after n_iter = 10 iteration steps. The surface is characterized by a very regular parameter grid without any artefacts. The described approach, therefore, yields very promising results. Further improvements can be obtained by examining the investigated data set for regularities in the elevation and to divide it into levels in such a way that the neighbouring point cloud profiles have a similar course. However, these are future investigations; in the following paragraphs equidistant profiles are used to construct skinned surfaces.

Figure 16:

Skinned surface for the data set Beach based on section curves that are determined by involving the neighbouring curves (n_iter = 10).

5.4 Local parameterization and iterative reparameterization

The skinned surface presented in Figure 16 represents the local structures of the beach quite well. The point cloud is already approximated very satisfactorily along the section curves, but, naturally, there are systematic deviations between the section curves. Furthermore, a large number of control points is required (the example surface in Figure 16 has 319 control points along the section curves). Besides, the parameterization reflects the local structures only in direction of the section curves, but not across them.

Due to these reasons, two different skinned surfaces are constructed: one in u-direction of the surface to be estimated and one in v-direction (S _u (u, v) and S _v (u, v), respectively). By projecting the point cloud onto the former one, the u-parameter is derived, whereas the latter one serves as base surface for parameterization in v-direction. Having derived the surface parameters, the best-fitting TP B-spline surface – the LP B-spline surface – can be estimated as described in Section 2.4 (cf. also Figure 10).

In analogy to Figures 5 and 17 (left) presents the data set Beach coloured according to the Euclidean residuals resulting from the approximation by means of an LP B-spline surface with 15 section curves. Compared to the results achieved by classical TP B-spline surfaces, a huge improvement can be observed: Rather than having maximal deviations of up to 10 m, they can be considerably reduced to a maximal magnitude of a decimeter. Compared to the measurement noise, these deviations are still loo large. However, the corresponding residuals (coloured in yellow in Figure 17 (left)) accumulate along lines that run parallel to the section curves and, thus, clearly indicate systematic deviations. The origin of these systematic deviations can be explained in particular by the most conspicuous yellow line in Figure 17 (left) that runs across the local maxima in the middle of the point cloud. The two adjacent parallel section curves both have a lower elevation than these local maxima, so that the maxima are “cut off” when constructing the control curves. Similar reasons cause the other stripes of large residuals. The location of the section curves, therefore, has a major influence on the approximation quality. This is also emphasized by Figure 17 (right), presenting the residuals resulting from a parameterization based on only seven section curves. Obviously, these seven curves catch the local extrema better than the 15 curves. As a consequence, the yellow lines are considerably less pronounced.

Figure 17:

Point cloud Beach coloured according to the Euclidean residuals [m] resulting from the approximation by means of LP B-spline surfaces with different numbers of section curves.

These results reveal that, rather than arbitrarily placing the section curves, care should be taken that they run through the local extrema of the point cloud. Yet, these investigations clearly are beyond the scope of this paper and will be postponed to the future.

Alternatively, the application of iterative reparameterization [33] with the skinned surfaces being only the initial solution for parameterization (cf. Figure 10) already considerably improves the approximation quality as demonstrated in Figure 18, presenting the residuals resulting from iterative reparameterization with an initial base surface constructed from 15 section curves.

Figure 18:

Point cloud Beach coloured according to the Euclidean residuals [m] resulting from iterative reparameterization with skinned surfaces serving as initial base surfaces.

Due to the large systematic residuals near the local extrema after the first iteration step, there are still accumulations of larger residuals in these regions after the final iteration step. However, they are not longer stripy, and they have been considerably reduced during reparameterization. Further reduction, if not even elimination, should be achieved by the targeted placement of the section curves to be investigated in future. Apart from these remaining slightly systematic deviations, Figure 18 indicates a random behaviour of the residuals within a magnitude that corresponds to the measurement accuracy of a long-range scanner.

The corresponding B-spline surface (see Figure 19) still has the compact form of a TP B-spline surface, is able to represent a variety of local structures and yet has a comparatively small number of control points (319 control points have been used for the skinned surfaces in each parameter direction, whereas for the final LP B-spline surface only 40 control points in each parameter direction have been used). Thus, this surface is perfectly suitable to serve as a basis for subsequent analysis steps.

Figure 19:

Approximating LP B-spline surface after iterative reparameterization.

6.1 Quality measures used for the analysis

The analysis of the approximation quality is based on three different quality measures:

1. The sum of squared residuals

   (11) Ω = ε T ε

   is minimized during the approximation procedure and evaluates the surface’s goodness-of-fit. In case Ω is compared for different B-spline surfaces with the same number of control points which are estimated by means of the same number of observations, the sum of squared residuals can be used in the sense of information criteria to choose the optimal model (see e.g. [38, 39]). Hence, Ω is used in Section 6.3 to compare the goodness-of-fit of TP and LP B-spline surfaces, with a smaller sum of squared residuals indicating the better model.

2. In order to derive information about systematic deviations between observations and estimated surface, the residuals’ probability distribution is analysed in Section 6.4. For this purpose, the residuals’ coordinate-wise histograms are regarded. In addition to a visual comparison, the first four statistical moments are computed and compared. For all coordinate directions, the histograms’ expectation value EV can be used to detect a biased estimation. Apart from that, different interpretations for the histograms of v_x,y and of v _z are necessary, as the estimated surface is oriented approximately parallel to the x-/y-plane for all investigated data sets. Thus, the allocated surface parameters have a direct influence on the x- and y-coordinate and, therefore, on the distribution of the corresponding residuals.

   For the histograms of v _z , a comparison of the standard deviations STD with the measurement uncertainty allows to evaluate the goodness-of-fit. Assuming normally distributed measurement noise, the histograms of v _z should be close to a normal distribution with a standard deviation that corresponds to the uncertainty of the laser scanner. Skewness SN and kurtosis KS indicate deviations from the normal distribution.

   Contrarily, a good parameterization decreases the residuals in x- and y-direction. Therefore, the more residuals v _x and v _y are close to zero, the better the parameterization and, thus, the approximation quality. Hence, for the residuals v _x and v _y , histograms with small standard deviations (even smaller than the measurement uncertainty) and a large kurtosis indicate a good approximation. Since systematics should be avoided also in these two coordinate directions, a skewness close to 0 is still aimed for.

3. Nominal surfaces are available for the data sets Specimen _u and Specimen _n . Thus, the estimated observations’ Euclidean distances d _i with respect to these nominal surfaces are regarded in Section 6.5. For a compact analysis, these distances are summarized by the mean value d ̄ as well as the maximal value d_max.

6.2 Settings for the estimation of B-spline surfaces

When estimating B-spline surfaces, different parameter groups have to be defined (cf. Section 2.4). For the investigations conducted in the following subsections, solely cubic B-spline surfaces (p = q = 3) are used.

In order to maintain the comparability between TP and LP B-spline surfaces, the same number of control points is used for both types of surfaces when approximating the same data set. Due to the loops and artefacts that emerge when using too large a number of control points for the TP B-spline surfaces, it may be that the chosen number of control points is not the optimal one for the LP B-spline surfaces. The choice of the optimal number of control points for LP B-spline surfaces, however, is not part of this publication, but will be investigated in future.

When using LP B-spline surfaces, the number of section curves in both parameter directions n_s,u and n_s,v, respectively, emerge as additional parameters to be chosen. In order to investigate these parameters’ influence on the approximation quality, different values are arbitrarily chosen for each estimated LP B-spline surface. The optimization of this parameter, however, is also part of future investigations.

Finally, a simplified stochastic model assuming independent and identically distributed observations is used. A complete stochastic model will be incorporated in future.

6.3 Analysis of the sum of squared residuals

Table 3 lists the sum of squared residuals for all data sets resulting from the modelling by means of TP B-spline surfaces and LP B-spline surfaces with varying number of section curves.

For all four data sets a clear improvement of the approximation quality in terms of the goodness-of-fit due to the local parameterization can be observed. Starting with the results of the non-noisy data set Specimen _u , the sum of squared residuals can be reduced by a factor of four when using a local parameterization that is based on n_s,u = 7 and n_s,v = 12 section curves. This result can be further improved by increasing the number of section curves. Although non of the investigated B-spline surfaces yields the perfect result of Ω = 0, which is achievable when using the nominal surface parameters, the results achieved by local parameterization come very close to it, even without optimizing the number, position and orientation of the section curves. The results of the noisy data set Specimen _n , which are slightly larger due to the superimposed noise, fully confirm the findings.

The reduction of the sum of squared residuals is most pronounced for the data set Beach, for which a reduction by more than a factor of 1000 can be observed. The improvement achieved for the data set Soil is much smaller, but nevertheless clearly noticeable. The large differences in the improvement of the sum of squared residuals for the data set Beach compared to the other data sets can be explained by the fact that the height variations, which limit the approximation by means of TP B-spline surfaces, are considerably more pronounced in the data set Beach (in the range of several decimeters to meters) than in the other three data sets (in the range of some millimeters to centimeters).

Increasing the number of section curves results in an improvement of the approximation quality in most cases. The only exception can be observed for the data set Beach, for which the use of n_s,u = n_s,v = 30 section curves slightly increases Ω compared to the use of n_s,u = n_s,v = 15 section curves. As already indicated by Figure 17, a disadvantegous placement of the section curves may cause this deterioration.

6.4 Analysis of the residuals’ histograms

In Figure 18, the Euclidean residuals’ spatial distribution has already been presented for the data set Beach, indicating a clear reduction of systematics compared to the residuals of the TP representation shown in Figure 5. These representations of the residuals are complemented by Figure 20, presenting the corresponding coordinate wise residuals’ histograms. In Figure 20 (left), the residuals in x-direction can be seen, which are also representative for the residuals in y-direction, with x and y being the main directions of the estimated surface. The histograms of the residuals in z-direction, being perpendicular to the estimated surface’s main directions, are shown in Figure 20 (right). Only a small improvement due to the local parameterization can be observed for the latter: The residuals’ histogram belonging to the LP B-spline surface (red histogram) is slightly narrower and has slightly more residuals close to zero than the blue histogram of the TP B-spline surface. In contrast, an enormous effect can be observed in the histogram of the residuals v _x : The histogram belonging to the TP B-spline surface is unexpectedly flat with a large number of residuals that have a magnitude of several centimeters and even decimeters (for reasons of clarity, the residuals’ range is restricted to [−0.2, 0.2] in Figure 20 (left), although there are residuals that exceed this order of magnitude). The shape of the histogram clearly indicates an unsatisfactory parameterization (cf. explanations in Section 6.1). Due to the local parameterization, the histogram is considerably narrowed, resulting in a significantly better distribution of the residuals.

Figure 20:

Histograms of the residuals resulting from the approximation of the data set Beach by means of a TP B-spline surface (30 ⋅ 30 control points) and by means of an LP B-spline surface (7 section curves, 30 ⋅ 30 control points).

The effect is this pronounced only for the data set Beach. It is closely related to the excessively large sum of squared residuals in Table 3 and, again, can be explained by the large height variations within the point cloud, which are far from being well represented by the initial base surface used for the estimation of a TP B-spline surface.

The characterizing parameters (expectation value EV, standard deviation STD, skewness SN, kurtosis KS) of the respective histograms are summarized in Table 4 for the same approximating B-spline surfaces that have been investigated in Table 3.

The analysis of the expectation vector reveals a minimally biased estimation for the TB B-spline surface, whereas the histograms belonging to the LP B-spline surfaces are perfectly centred around 0. The standard deviation in z-direction can be reduced by local parameterization from 2.7 cm to 1.5 cm, which perfectly well coincides with the measurement uncertainty of a long-range laser scanner. Skewness and kurtosis of the histograms in z-direction indicate a clear deviation of the residuals’ distribution from a normal distribution when using TP B-spline surfaces. A convergence towards a normal distribution can be observed due to the local parameterization. However, skewness and kurtosis both still deviate from the target values of a normal distribution and, hence, indicate remaining systematics in the residuals. They are likely to be further reduced by optimizing the position, orientation and number of section curves.

The results of the histograms in x- and y-direction are similar promising. As already expected from Figure 20, the biggest improvement due to local parameterization can be achieved with respect to the residuals’ standard deviations: The standard deviations in x- (σ _x ≈ 30 cm) and y-direction (σ _y ≈ 55 cm) are disproportionately large for the TP B-spline surface. After local parameterization, they are with < 6 mm clearly smaller than the standard deviation in z-direction and, thus, indicate a considerably improved parameterization (cf. explanations in Section 6.1).

The histograms’ kurtosis is clearly increased by local parameterization for the residuals v_x,y, resulting in leptocurtic histograms. This effect is most pronounced for the residuals in y-direction. Hence, the residuals accumulate near zero as requested (cf. explanations in Section 6.1). Solely the skewness does not show the desired behaviour for all histograms. For two of the LP B-spline surfaces, an increase of the skewness is observable compared to the results obtained by the TP B-spline surface. Obviously, systematics are introduced into the residuals by these two LP B-spline surfaces. However, the fact that the third LP B-spline surface reduces the skewness, reinforces the already gained impression that the position and number of the section curves must be chosen very carefully, especially for such challenging data sets as the data set Beach.

For the data set Beach, the unknown and, therefore, unmodelled stochastic properties of the laser scanner used as well as the scanning configuration, in particular the angle of incidence, are likely to influence the parameters presented in table 4. Therefore, the results are complemented by Table 5, presenting the corresponding results of the data set Specimen _n , for which these influences do not occur.

For this data set, the histogram’s parameter in z-direction are already close to the target values when using a TP B-spline surface: The standard deviation and the kurtosis are only slightly too large, whereas the skewness indicates some remaining systematics. Nevertheless, these results can be further improved by local parameterization. All three LP B-spline surfaces produce histograms that almost perfectly reproduce the properties of the simulated measurement noise. The parameters describing the histograms of v_x,y also confirm the findings derived by means of Table 4.

6.5 Comparison with nominal surfaces

Nominal surfaces, which allow the evaluation of the estimations’ correctness, are available for the data sets Specimen _u and Specimen _n . As an example, Table 6 lists the mean and maximal Euclidean deviations with respect to the nominal surface for data set Specimen _u . Particularly striking is the large maximal deviation that results from the approximation by means of the TP B-spline surface: With almost 1.3 cm this value is more than four times larger than the corresponding values resulting from the approximation by means of LP B-spline surfaces. In order to investigate the origin of the large discrepancies, Figure 21 presents the estimated observations of Specimen _u coloured according to the Euclidean distances [mm] with respect to the nominal surface for TP_10⋅10 and LP_9,17,10⋅10. When comparing the two figures, two things are particularly striking. Firstly, comparatively large deviations ( > 2.5 mm ) with respect to the nominal surface occur in the boundary regions when using the TP B-spline surface. These large deviations are caused by the boundary constraints that are used to prevent the degeneration of the TP B-spline surface during the iterative reparameterization (cf. [33]). As the LP B-spline surfaces are much more stable during reparameterization, no constraints are required and, hence, no accumulations of large deviations occur in the boundary regions. These accumulations are the reason for the large maximal and mean deviation produced by the TP B-spline surface in Table 6. In the inner part of the surface, however, the deviations for TP and LP B-spline surfaces attain values of similar magnitudes. Obviously, the data set Specimen _u is due to its small height changes considerably more good-natured with respect to the approximation by means of TP B-spline surfaces than e.g. the data set Beach. The second striking feature in Figure 21 is the behaviour of individual estimated observations that are obviously shifted out of the actual sampling pattern when using the TP B-spline surface, resulting in small white holes in the point cloud presented in Figure 21 (left). This effect does not occur with the LP surface, rather the sampling pattern is perfectly preserved during parameterization.

Figure 21:

Estimated observations of Specimen _u coloured according to the Euclidean distances ([mm]) with respect to the nominal surface. In order to retain comparability and in order to reveal patterns, the colorbars are restricted to a maximum deviation of 2.5 mm.

7.1 Summary

In this contribution, an approach was developed to improve the approximation quality in terms of goodness-of-fit, correctness and stability of classical TP B-spline surfaces. As the surface parameters contain a huge amount of information regarding the surface to be estimated, the basic idea of the developed strategy is to improve the initial parameterization, which is then further improved by iterative reparameterization. The improvement of the initial parameterization is achieved by using two skinned B-spline surfaces (one in each parameter direction of the surface to be estimated) as initial base surfaces. As the skinned B-spline surfaces are constructed by means of (almost) independently determined B-spline curves, which are able to represent local undulations, the initial parameterization already reproduces local structures that are contained within the point cloud to be approximated. Here, it turned out that a prerequisite for the successful use of skinned surfaces as base surfaces is the involvement of the neighbouring curves in the estimation of the individual section curves as, otherwise, discontinuities emerge in the parameter space of the skinned surface.

The developed LP B-spline surfaces were then used to model four data sets (simulated and measured) with different characteristics, and the results were compared to nominal surfaces as well as to the results achieved by means of classical TP B-spline surfaces. For all investigated data sets, an improvement due to the local parameterization in terms of goodness-of-fit, correctness and the resulting residuals’ distribution can be observed. For good-natured data sets with rather small height changes, the improvement of the considered quality measures is noticeable, but not overly large. Not without reason are the classical TP B-spline surfaces already used in many applications. However, in addition to the small quantitative improvement due to the local parameterization, the process of iterative reparameterization is considerably stabilized when starting with the improved base surfaces. As a consequence, no constraints have to be introduced into the adjustment to prevent the surface’s degeneration. For more challenging data sets (in this contribution the data set Beach), which are characterised by large height changes, the improvement is tremendous: While these data sets cannot at all be satisfactorily approximated by means of TP B-spline surfaces, a very satisfying approximation is yielded by the developed LP B-spline surfaces.

Due to the storage of the relevant information regarding the surface to be estimated in the surface parameters, LP B-spline surfaces get along with a comparatively small number of control points. Solely the intermediate results, the skinned base surfaces, require a relatively large number of control points. Furthermore, unlike with refinement strategies like LR or HB-splines, the compact form of TP B-spline surfaces is maintained. Both properties enable a simple use of LP B-spline surfaces in further processing steps like deformation analyses.

7.2 Future investigations

The approximation results yielded by the developed LP B-spline surfaces are already very promising. At the same time, however, they indicate potential for improvement. The results reveal that the position of the section curves in particular has a large influence on the quality of the point cloud modelling. Closely linked to the section curves’ position is their orientation and their number. Hence, strategies to automatically determine these parameters are going to be developed in future.

Furthermore, the data sets investigated are currently limited in one respect, as they are restricted to have a rectangular boundary. The extension of the developed approach to point clouds with arbitrary boundaries also requires further investigations.

Another important aspect to be investigated in future is the differentiation between outliers and locally differing surface patters, which becomes particularly important as soon as accumulating outliers are contained in the point clouds to be modelled.

Finally, the comparisons conducted in this paper are exclusively made with TP B-spline surfaces. In the future, these comparisons will be supplemented by comparisons with refinement strategies like LR B-splines and HB-splines. The comparison will be made both in terms of the approximation quality and with regard to computational demands.