Laser Scanner–Based Deformation Analysis Using Approximating B-Spline Surfaces

Abstract

Due to the increased use of areal measurement techniques, such as laser scanning in geodetic monitoring tasks, areal analysis strategies have considerably gained in importance over the last decade. Although a variety of approaches that quasi-continuously model deformations are already proposed in the literature, there are still a multitude of challenges to solve. One of the major interests of engineering geodesy within monitoring tasks is the detection of absolute distortions with respect to a stable reference frame. Determining distortions and simultaneously establishing the joint geodetic datum can be realised by modelling the differences between point clouds acquired in different measuring epochs by means of a rigid body movement that is superimposed by distortions. In a previous study, we discussed the possibility of estimating these rigid body movements from the control points of B-spline surfaces approximating the acquired point clouds. Alternatively, we focus on estimating them by means of constructed points on B-spline surfaces in this study. This strategy has the advantage of larger redundancy compared to the control point–based strategy. Furthermore, the strategy introduced allows for the detection of rigid body movements between point clouds of different epochs and for the simultaneous localisation of areas in which the rigid body movement is superimposed by distortions. The developed approach is based on B-spline models of epoch-wise acquired point clouds, the surface parameters of which define point correspondences on different B-spline surfaces. Using these point correspondences, a RANSAC-approach is used to robustly estimate the parameters of the rigid body movement. The resulting consensus set initially defines the non-distorted areas of the object under investigation, which are extended and statistically verified in a second step. The developed approach is applied to simulated data sets, revealing that distorted areas can be reliably detected and that the parameters of the rigid body movement can be precisely and accurately determined by means of the strategy.

1. Introduction

Geodetic deformation monitoring deals with the measurement, modelling and evaluation of geometric object changes caused by external influences, such as temperature changes, load effects or changes in the groundwater level [1]. Classically, point-based strategies are used to monitor the behaviour of objects over time: By repeatedly measuring signalised characteristic points of the object under investigation, displacement vectors can be derived, representing the object’s deformation [2]. Due to its long tradition in engineering geodesy, sophisticated strategies to analyse point-wise deformation measurements are established. However, point-based deformation analysis also has a variety of drawbacks: In order to represent the behaviour of the entire object of interest by means of single points, characteristic points have to be carefully chosen, often by including prior knowledge of the expected deformation. Regardless of a successful discretisation process, the derived deformation measures are always sparse. Furthermore, the repeated measurement of the same characteristic points requires a signalisation of these points and, thus, access to the object under investigation. Finally, point-based approaches can be very time and labour intensive, especially when large objects are monitored [3,4,5].

The terrestrial laser scanner overcomes these drawbacks of point-based measuring strategies: a laser scanner remotely samples even large objects with a high spatial resolution in a very short time and directly provides a three-dimensional and quasi-continuous description of the object under investigation [3,6,7]. However, despite the metrological strengths of laser scanning, the analysis of the resulting point clouds—especially with regard to a deformation analysis—poses a variety of challenges. Among other things, the measuring principle of laser scanners does not usually allow for the direct reproduction of points in subsequently acquired point clouds [6] and, thus, the direct computation of displacement vectors between different point clouds is not expedient. Furthermore, the single point precision within a laser scanning point cloud may be considerably lower than, for example, the precision of tacheometrically measured signalised points [8]. Nevertheless, numerous strategies to compare point clouds and, thus, to perform a point cloud–based deformation analysis are proposed in the literature. These strategies can be categorised into three classes [6]:

• When performing a point-to-point (P2P) comparison, it is assumed that point correspondences in different point clouds either exist or can be constructed. Based on these identical points, displacement vectors can be determined. In some situations, the required point correspondences can be established by an appropriate measurement setup (e.g., [9,10]). Alternatively, they can either be constructed by means of local models of at least one of the point clouds (e.g., [11,12]) or by means of feature descriptors (e.g., [13,14]).
• When implementing a point-to-surface (P2S) comparison, one of the point clouds—usually the point cloud of the first measuring epoch—is approximated by a reference surface. In order to analyse point clouds acquired in subsequent epochs for deformations, the distance with respect to the reference surface is determined for each point of these point clouds. Depending on the type of reference surface chosen, mesh-based approaches (e.g., [15,16]) and approaches based on analytical surfaces (e.g., [12,17]) can be distinguished.
• When using a surface-to-surface (S2S) comparison, the entirety of the point clouds is modelled either by means of meshes or by means of analytical surfaces. The examination for deformations between these models can be performed in two different ways. The first form of implementation uses the models of the point clouds to generate identical points in the different epochs, allowing for the computation of displacement vectors (e.g., [18,19]). The second form of implementation is based solely on analytical models of the point clouds and compares the estimated parameters of the respective models (e.g., [10,20]).

Strategies using models of point clouds take advantage of the high redundancy in the point clouds in order to overcome the challenge of the comparatively low single point precision within laser scanning point clouds: The theoretical precision of a surface that is determined on the basis of these points is significantly higher than that of the single points [8]. However, two more challenges can be defined when using a model-based analysis approach [21]: Firstly, there is still a lack of interpretable deformation measures between two surfaces. Although displacement vectors fulfil this property, a prerequisite for their use is the correct definition of point correspondences on the surfaces. Secondly, in only a few of the approaches presented in the literature, the deformation’s significance is confirmed by means of statistical tests.

Freeform surfaces, such as B-splines, allow for the modelling of different kinds of objects with various geometric properties. Being parametric surfaces, they directly provide the possibility to define point correspondences on different B-spline surfaces, assuming that the surface parameters are consistently chosen. Moreover, their property of locality allows for the simultaneous modelling of distorted regions of an object and its non-distorted parts by means of one single B-spline surface.

The aim of the present publication is the development of a B-spline-based strategy for performing a point cloud–based deformation analysis, using these beneficial properties of B-splines. Unlike most of the strategies available in the literature, our approach simultaneously determines rigid body movements and distortions of the object under investigation. Being designed as an S2S comparison, all acquired point clouds are modelled by means of consistently parameterised B-spline surfaces, the surface parameters of which are used to define point correspondences. Based on these point correspondences, rigid body movements between the point clouds are estimated, using a robust random sample consensus (RANSAC) approach. The resulting consensus set defines an initial region that is not subjected to additional distortions. In the final step, these stable regions are extended and statistically verified, providing the basis for the computation of displacement vectors in the distorted areas. Both the estimated rigid body movements and the identified distorted regions are confirmed by means of statistical hypothesis testing. Hence, the developed approach combines the strength of areal analysis strategies with the sophisticated tools of the classical point-based deformation analysis.

2. Materials and Methods

2.1. Estimation of Rigid Body Movements and Simultaneous Localisation of Non-Distorted Regions Using Sampled Approximating B-Spline Surfaces

Commonly, the term ‘deformation’ summarises two types of geometric object changes: an object underlying rigid body movements is subject exclusively to rotations and translations, whereas an object underlying distortions changes in shape [1]. The aim of the strategy introduced in the following subsections is to robustly estimate rigid body movements of an object and to simultaneously localise the object’s distorted parts, providing the basis for the determination of displacement vectors.

2.1.1. B-Spline Surfaces for Deformation Analysis

The developed approach is based on approximating B-spline surfaces. The mathematical definition of a B-spline surface is given by the following [22]:

$$\widehat{\mathbf{S}}(u,v)=\left[\begin{array}{c}{\widehat{S}}_x(u,v)\\ {\widehat{S}}_y(u,v)\\ {\widehat{S}}_z(u,v)\end{array}\right]=\sum _{i=0}^n\sum _{j=0}^m{N}_{i,p}\left(u\right)N_{j,q}\left(v\right){\mathbf{P}}_{ij},\begin{array}{c}\hfill 0\le u,v\le 1\end{array}.$$

(1)

According to Equation (1), a three-dimensional surface point $\widehat{\mathbf{S}}(u,v)$ in Cartesian coordinates is computed as the weighted average of $(n+1)·(m+1)$ control points ${\mathbf{P}}_{ij}$. The B-spline basis functions $N_{i,p}\left(u\right)$ and $N_{i,q}\left(v\right)$ with degree p and q are the corresponding weights and can be computed recursively (cf. references [23,24]). The basis functions are functions of the surface parameters u and v, locating a surface point on the surface and defining the parameter space. Two knot vectors $\mathbf{U}=[u_0,\dots ,u_r]$ and $\mathbf{V}=[v_0,\dots ,v_s]$ split the B-splines’ domain into knot spans, leading to the B-splines’ property of ‘locality’, meaning that the shifting of a single control point changes the surface only locally.

B-spline surfaces allow for the flexible modelling of a variety of objects under investigation. Hence, they are used in several engineering geodetic applications, either for geometric state descriptions (cf. references [25,26,27]) or in the context of monitoring tasks ([28,29]). In addition to their flexibility, B-spline surfaces have two other properties that qualify them for use within an areal deformation analysis:

• Firstly, B-spline surfaces are invariant under common geometric transformations, such as the similarity transform: the application of a transformation to the control points ${\mathbf{P}}_{ij}$ of a B-spline surface results in the same surface as the application of this transformation to the surface itself [22]. As proven in [20], the estimated control points of different B-spline surfaces can be used as corresponding points to estimate the parameters of a similarity transform, provided that a joint parameterisation is used when determining the best-fitting B-spline surfaces. However, not only can the control points be used for this purpose, but also surface points $\widehat{\mathbf{S}}(u,v)$, with the surface parameters u and v identifying the point correspondences. This property is demonstrated in Figure 1, presenting two B-spline surfaces that differ from each other only by rotations and translations (left surface: epoch 1, right surface: epoch 2). As can be seen, parameter grids (green quadrangles) of the surface of epoch 1 that are individually transformed by means of the same transformation parameters as the surface are congruent with the respective grid quadrangles of the transformed surface.
• Secondly, the locality of the B-spline surfaces directly allows for the modelling and for the localisation of an object’s distortions. For example, when the rigid body movement of a B-spline surface is superimposed by distortions, non-distorted areas can be identified by means of the surface parameter grid as demonstrated in Figure 2. The figure presents two B-spline surfaces that differ from each other by a rigid body movement and a superimposed local distortion (left surface: epoch 1, right surface: epoch 2). The transformation of the four parameter quadrangles (green quadrangles) of the surface of epoch 1, using the same transformation parameters as those for the rigid body movement of the surface, highlights that this local distortion only locally influences the parameter lines: only the red quadrangle in epoch 2—being the only one of the four example quadrangles located in the distorted area—is not congruent with the respective grid quadrangles of the surface of epoch 2.

These two properties of B-spline surfaces form the cornerstones of the developed approach that is summarised in the flow chart in Figure 3.

For the sake of simplicity, only two measuring epochs are considered; the extension by further epochs is straightforward.

Starting points are point clouds of an object acquired in two different measuring epochs ($P{C}^{\left(1\right)}$ and $P{C}^{\left(2\right)}$), along with their variance covariance matrices (VCM) (${\mathsf{\Sigma }}_{P{C}^{\left(1\right)}}$ and ${\mathsf{\Sigma }}_{P{C}^{\left(2\right)}}$). Between these two epochs, the object may undergo rigid body movements, which additionally may be superimposed by distortions. It is assumed, however, that only parts of the object are subject to distortions. The resulting non-distorted regions of the object are necessary for successful application of the developed approach.

Both point clouds are initially approximated by means of best-fitting B-spline surfaces, taking into account the stochastic information of the available point clouds (Section 2.1.2). Due to the use of joint parameterisation during this step, afterwards, the surface parameters can be used to defined point correspondences on the estimated B-spline surfaces. These corresponding points, resulting from discretisations of the B-spline surfaces (Section 2.1.3), are then used to robustly estimate rigid body movements and to initially define non-distorted regions of the object (Section 2.1.4 and Section 2.1.5). Finally, the non-distorted regions are extended and their stability as well as their extents are statistically confirmed (Section 2.1.6).

2.1.2. Point Cloud Modelling by Means of Best-Fitting B-Spline Surfaces

When using B-spline surfaces (cf. Equation (1)) to model a point cloud $P{C}^{\left(i_e\right)}$ ($i_e=1,2$: indicating the measuring epoch), the $n_l^{\left(i_e\right)}$ individual points of the point cloud are the observations ${\mathbf{S}}_k^{\left(i_e\right)}(u,v)$ ($k=1,\dots ,n_l^{\left(i_e\right)}$) that are used to determine the best-fitting B-spline surface ${\widehat{\mathbf{S}}}_k^{\left(i_e\right)}(u,v)$:

$$\begin{array}{cc}\hfill {\widehat{\mathbf{S}}}_k^{\left(i_e\right)}(u,v)& ={\mathbf{S}}_k^{\left(i_e\right)}(u,v)+{\epsilon }_k^{\left(i_e\right)}(u,v),\hfill \end{array}$$

(2)

with the three-dimensional residual vectors ${\epsilon }_k^{\left(i_e\right)}(u,v)$. As in Equation (1), each observation ${\mathbf{S}}_k^{\left(i_e\right)}(u,v)$ is a three-dimensional point expressed in Cartesian coordinates.

Usually, only the position of the control points ${\mathbf{P}}_{ij}$ (cf. Equation (1), for reasons of simplicity, the epochal affiliations indicated by the superscript $\left(i_e\right)$ are not carried along) is estimated in a linear Gauß Markov model, according to the following [31]:

$$\begin{array}{cc}\hfill {\widehat{\vartheta }}_{\mathbf{P}}& ={({\mathbf{A}}^T·{\mathbf{Q}}_{\mathbf{ll}}^{-1}·\mathbf{A})}^{-1}{\mathbf{A}}^T·{\mathbf{Q}}_{\mathbf{ll}}^{-1}·\mathbf{l}.\hfill \end{array}$$

(3)

In Equation (3), the vector of unknowns ${\widehat{\vartheta }}_{\mathbf{P}}$ summarises the estimated positions of all $(n+1)·(m+1)$ control points and is ordered coordinate-wise. The observation vector $\mathbf{l}$ is structured accordingly and contains the observed noisy surface points ${\mathbf{S}}_k(u,v)$. The stochastic behaviour of the observations is characterised by the corresponding cofactor matrix ${\mathbf{Q}}_{\mathbf{ll}}$, the inverse of which is used as a weighting matrix in the estimation and which differs from the corresponding VCM by the a priori variance factor ${\sigma }_0^2$ [1]:

$$\begin{array}{cc}\hfill {\mathsf{\Sigma }}_{\mathbf{ll}}& ={\sigma }_0^2·{\mathbf{Q}}_{\mathbf{ll}}.\hfill \end{array}$$

(4)

The design matrix $\mathbf{A}$, describing the functional relationship between observations and unknowns, is determined by means of the B-spline basis functions [31]:

$$\mathbf{A}=\left[\begin{array}{ccc}{N}_{0,p}\left(u_1\right)N_{0,q}\left(v_1\right)& \dots & N_{n,p}\left(u_1\right)N_{m,q}\left(v_1\right)\\ ⋮& \ddots & ⋮\\ N_{0,p}\left(u_{n_l}\right)N_{0,q}\left(v_{n_l}\right)& \dots & N_{n,p}\left(u_{n_l}\right)N_{m,q}\left(v_{n_l}\right)\end{array}\right]\otimes {\mathbf{I}}_{3\times 3},$$

(5)

with ${\mathbf{I}}_{3x3}$ being the 3 × 3 identity matrix. The vector of estimated control points ${\widehat{\vartheta }}_{\mathbf{P}}$ is completed by the corresponding VCM ${\mathsf{\Sigma }}_{\widehat{\vartheta }\widehat{\vartheta }}$, containing the stochastic information of the estimated control points [1]:

$$\begin{array}{cc}\hfill {\mathsf{\Sigma }}_{\widehat{\vartheta }\widehat{\vartheta }}& ={\widehat{\sigma }}_0^2·{({\mathbf{A}}^T·{\mathbf{Q}}_{\mathbf{ll}}^{-1}·\mathbf{A})}^{-1},\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}{\widehat{\sigma }}_0^2=\frac{{\widehat{\epsilon }}^T·{\mathbf{Q}}_{\mathbf{ll}}^{-1}·\widehat{\epsilon }}{3(n_l-(n+1)·(m+1))}\hfill \end{array}$$

(6)

and with the vector of estimated residuals $\widehat{\epsilon }$ being computed according to the following:

$$\begin{array}{cc}\hfill \widehat{\epsilon }& =\widehat{\mathbf{l}}-\mathbf{l}=\mathbf{A}{\widehat{\vartheta }}_{\mathbf{P}}-\mathbf{l}.\hfill \end{array}$$

(7)

A geometrical interpretation of the estimation $\widehat{\mathbf{l}}=\mathbf{A}{\widehat{\vartheta }}_{\mathbf{P}}$ in Equation (7) is that the original vector of observations $\mathbf{l}$ is orthogonally projected to the $3(n+1)·(m+1)$-dimensional column space of the design matrix $\mathbf{A}$ [32]. Thus, the vector of estimated observations $\widehat{\mathbf{l}}$ can be computed as a linear combination of the vector of estimated control points ${\widehat{\vartheta }}_{\mathbf{P}}$. Due to $n_l>(n+1)·(m+1)$, linear dependencies between the estimated observations emerge. The corresponding VCM of the vector of estimated observations

$$\begin{array}{cc}\hfill {\mathsf{\Sigma }}_{\widehat{\mathbf{l}}\widehat{\mathbf{l}}}& =\mathbf{A}·{\mathsf{\Sigma }}_{\widehat{\vartheta }\widehat{\vartheta }}·{\mathbf{A}}^T,\hfill \end{array}$$

(8)

thus, has a maximum rank of $rk\left({\mathsf{\Sigma }}_{\widehat{\mathbf{l}}\widehat{\mathbf{l}}}\right)=3(n+1)·(m+1)$ and—assuming an overdetermined estimation problem—is singular.

The strategy described above requires appropriately determining the remaining parameter groups of the best-fitting B-spline surface separately from the adjustment of the optimal control points: The B-spline’s degrees p and q are usually specified a priori. Using cubic B-splines with $p=3$ and $q=3$ is a generally accepted choice [22]. Strategies to determine the knot vectors are proposed in, for example, [22,31,33]. The optimal number of control points to be estimated can be interpreted as a model selection problem and can either be solved by classical model selection strategies or by structural risk minimisation (cf. [34,35]). The latter strategy has the advantage that the degrees p and q can also be taken into account during model selection, rather than choosing them arbitrarily [36]. Finally, appropriate surface parameters u and v have to be allocated to the observations (cf. e.g., [37] for an iterative parameterisation approach for laser scanning point clouds). When B-spline surfaces are used as a basis for deformation analysis, however, an independent parameterisation of different point clouds is not expedient; rather, a joint parameterisation has to be implemented as demonstrated in [35]. As simulated data sets are investigated in this initial study (cf. Section 2.2), nominal values for degrees, knot vectors and surface parameters are available and, therefore, the appropriate determination of these parameter groups can be left aside during this initial study.

2.1.3. Constructing Identical Points on Approximating B-Spline Surfaces

Having estimated the best-fitting B-spline surfaces, they are used to construct $n_{id}=n_{id,u}·n_{id,v}$ identical (id) points ${\mathbf{X}}_{k,l}^{\left(i_e\right)}$ in the two measuring epochs $i_e=1$ and $i_e=2$. Provided that the surfaces are based on a joint parameterisation, the surface parameters u and v, discretising the estimated surfaces (cf. Equations (1) and (2)), can be used for this purpose:

$$\begin{array}{cc}\hfill {\mathbf{X}}_{k,l}^{\left(i_e\right)}& ={\widehat{\mathbf{S}}}^{\left(i_e\right)}(u_k,v_l),\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h} k=1,...,n_{id,u},l=1,...,n_{id,v}\mathrm{a}\mathrm{n}\mathrm{d} i_e=1,2.\hfill \end{array}$$

(9)

Using the stochastic information of the estimated control points, the sampled surface points are completed by their VCMs ${\mathsf{\Sigma }}_{{\mathbf{X}}^{\left(i_e\right)},{\mathbf{X}}^{\left(i_e\right)}}$. Following the considerations in Section 2.1.2, linear dependencies emerge between the discretised points as soon as $n_{id}>(n+1)·(m+1)$, resulting in rank defects of the corresponding VCMs ${\mathsf{\Sigma }}_{{\mathbf{X}}^{\left(i_e\right)},{\mathbf{X}}^{\left(i_e\right)}}$. More precisely, a rank defect occurs in at least one parameter direction of the B-spline surface when the number of discretised surface points exceeds the number of estimated control points. The rank defect linked to the u-direction, therefore, can be computed according to the following [30]:

$${\delta }_u=\left\{\begin{array}{cc}0\hfill & \mathrm{if}{n}_{id,u}\le n+1\hfill \\ 3·(n_{id,u}-(n+1))(m+1)\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e},\hfill \end{array}\right.$$

(10)

and the rank defect linked to the v-direction accordingly arises to the following:

$${\delta }_v=\left\{\begin{array}{cc}0\hfill & \mathrm{if}{n}_{id,v}\le m+1\hfill \\ 3·(n_{id,v}-(m+1))(n+1)\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}.\hfill \end{array}\right.$$

(11)

The overall rank defect $\delta$ is the sum of these two parts [30]:

$$\begin{array}{cc}\hfill \delta & ={\delta }_u+{\delta }_v.\hfill \end{array}$$

(12)

This singularity in the VCMs of the corresponding points represents an important difference to deformation models in which the observed data points are directly incorporated.

2.1.4. Estimation of Rigid Body Movements

The estimation of rigid body movements of an object of interest is based on $n_{id}$ three-dimensional corresponding points ${\mathbf{X}}_{k,l}^{\left(1\right)}$ and ${\mathbf{X}}_{k,l}^{\left(2\right)}$ ($k=1,\dots ,n_{id,u}$, $l=1,\dots ,n_{id,v}$) in two measuring epochs, in this study on corresponding surface points (cf. Equation (9)).

Assuming that the object undergoes solely rigid body movements, all corresponding points—just like the entire surface—are subject to translations (summarised in the translation vector $\mathbf{t}$) and rotations (fully described by means of the rotation matrix $\mathbf{R}$, containing the rotation angles $\omega$, $\varphi$ and $\kappa$). Sometimes, a change of scale (indicated by the scale factor s) is also included, resulting in the mathematical formulation of a three-dimensional similarity transform [1]:

$$\begin{array}{cc}\hfill {\mathbf{X}}_{k,l}^{\left(2\right)}& =s·\mathbf{R}·{\mathbf{X}}_{k,l}^{\left(1\right)}+\mathbf{t},k=1,\dots ,n_{id,u},l=1,\dots ,n_{id,v}.\hfill \end{array}$$

(13)

When estimating the parameters of the similarity transform $\mathbf{t}$, $\mathbf{R}$ and s based on the corresponding surface points, it must be taken into account that both the points ${\mathbf{X}}_{k,l}^{\left(1\right)}$ available in the start system and the points ${\mathbf{X}}_{k,l}^{\left(2\right)}$ available in the target system are defined by estimated B-spline surfaces and, thus, are subject to uncertainties. In order to include the stochastic information of both point groups into the estimation of the transformation parameters, two strategies, delivering identical results, exist [38]: Either the transformation parameters are estimated in a Gauß–Helmert model, or the classical Gauß–Markov model, assuming the points in the start system to be deterministic, is extended. Because of the easier possibility to extend the approach to a robust one, the second variant is chosen here [30].

Assuming that the data sets under investigation are acquired by means of the same laser scanner and, hence, neglecting the scale factor s, the extended functional model to estimate the transformation parameters is then given by the following [38]:

$$\begin{array}{cc}\hfill {\widehat{\mathbf{X}}}_{k,l}^{\left(2\right)}& ={\mathbf{X}}_{k,l}^{\left(2\right)}+{\mathbf{e}}_{k,l}^{\left(2\right)}=\mathbf{R}·{\mathbf{X}}_{k,l}^{*\left(1\right)}+\mathbf{t}\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill {\widehat{\mathbf{X}}}_{k,l}^{\left(1\right)}& ={\mathbf{X}}_{k,l}^{\left(1\right)}+{\mathbf{e}}_{k,l}^{\left(1\right)}={\mathbf{X}}_{k,l}^{*\left(1\right)},\hfill \end{array}$$

(15)

with ${\mathbf{e}}_{k,l}^{\left(i_e\right)}$ ($i_e=1,2$) being residual vectors. The first part of the model (Equation (14)) describes the functional relationship between the identical points in the two epochs by means of the unknown transformation parameters. Equation (15) simultaneously introduces the coordinates of the identical points in the start system as additional observations and, thus, allows for the consideration of their stochastic information. The link between Equations (14) and (15) is given by the estimated point coordinates available in the start system, which are introduced in both equations and are marked with an asterisk. Neglecting inter-epochal correlations, the extended stochastic model is composed of the VCMs ${\mathsf{\Sigma }}_{{\mathbf{X}}^{\left(1\right)}{\mathbf{X}}^{\left(1\right)}}$ as well as ${\mathsf{\Sigma }}_{{\mathbf{X}}^{\left(2\right)}{\mathbf{X}}^{\left(2\right)}}$ of the constructed surface points, resulting in the following [38]:

$${\mathsf{\Sigma }}_{\mathbf{l}\mathbf{l}}=\left[\begin{array}{cc}{\mathsf{\Sigma }}_{{\mathbf{X}}^{\left(2\right)}{\mathbf{X}}^{\left(2\right)}}& \mathbf{0}\\ \mathbf{0}& {\mathsf{\Sigma }}_{{\mathbf{X}}^{\left(1\right)}{\mathbf{X}}^{\left(1\right)}}\end{array}\right]={\sigma }_0^2·{\mathbf{Q}}_{\mathbf{ll}}.$$

(16)

Based on the functional model in Equations (14) and (15) as well as on the stochastic model in Equation (16), the unknown transformation parameters can be estimated in a non-linear Gauß–Markov model (cf. [38] for more information).

2.1.5. Robust Estimation of Rigid Body Movements

The derivations in Section 2.1.4 are based on the assumption that the object under investigation is only subject to rigid body movements. However, rigid body movements are usually superimposed by additional distortions in reality. When taking into account points that are subject to both rigid body movements and distortions in the estimation of the rigid body movement, these points falsify the estimated transformation parameters. As a consequence, the determination of transformation parameters must necessarily be accompanied by an elimination of surface points in distorted areas. Otherwise, the required congruence of identical points is not given.

Considering points that are subject to both rigid body movements and distortion to be outliers in the context of the similarity transform is one possible strategy to cope with this challenge [30]. In order to accurately estimate rigid body movements in this situation, a robust strategy based on the random sample consensus (RANSAC)-algorithm [39] is proposed in this section. In general, RANSAC is an easy-to-implement method and is characterised by a high breakpoint (>50%). Furthermore, with respect to the presented problem, RANSAC allows for the robust estimation of rigid body movements while obtaining initial information about non-distorted regions.

The general idea of RANSAC is to initially determine the unknown parameters of a model describing $n_l$ data points by the minimum number $n_{sub}$ of data points that is needed for this determination [39]. These $n_{sub}$ data points are randomly drawn from the entire set of data points. After having computed the model parameters by means of these $n_{sub}$ points, the distances of the resulting preliminary model to the entirety of the data points are calculated. Points with distances that are within a predefined error tolerance $ϵ$ are considered to be consistent with the determined model parameters. These points form the consensus set ${\mathbb{S}}_{CS}$. The procedure is repeated iteratively until either the number of points in the consensus set equals a specified minimum number $n_{min}$ or a maximum number of iterations $i_{max}$ is performed. Finally, the points in the largest consensus set are used to estimate the optimal model parameters.

RANSAC can be directly applied to the robust estimation of rigid body movements that are superimposed by distortions (see Figure 4 for a schematic sketch of the procedure): the six parameters of a three-dimensional similarity transform with neglected scale factor (cf. Equation (14)) can principally be determined by means of two pairs of three-dimensional corresponding data points. However, the estimation of the transformation parameters described in Section 2.1.4 is a non-linear estimation problem and, thus, requires approximate values. As these approximate values are determined by means of quaternions (see, for example, [40]) in this study, a strategy that also yields the scale factor s, the number of point pairs required is $n_{sub}=3$. Hence, three pairs of sampled surface points are repeatedly and randomly drawn from the entirety of surface points and are afterwards used to determine approximate parameters ${\mathbf{t}}_0$, ${\mathbf{R}}_0$ and $s_0$ by means of a quaternion approach. (Strictly speaking, the modelled transformation is not a rigid body movement in this step. However, as the results are only approximate values for the final estimation of rigid body movements, this strategy is not to be regarded as critical.) Unlike the quaternion approach, the strategy presented in Section 2.1.4 allows for the accuracies of the corresponding points to be taken into account when estimating rigid body movements. For this reason, the approximate parameters are improved in a second step as described above, yielding optimal translations ${\widehat{\mathbf{t}}}_0$ and optimal rotations ${\widehat{\mathbf{R}}}_0$, describing the rigid body movements of the three chosen surface points. During this step, the scale factor is set to $\widehat{s}=1$ [30]. The entirety of computed surface points of the first measuring epoch ${\mathbf{X}}_{k,l}^{\left(1\right)}$ is afterwards transformed using the estimated transformation parameters as follows:

$$\begin{array}{cc}\hfill {\mathbf{X}}_{tr,k,l}^{\left(1\right)}& ={\widehat{\mathbf{R}}}_0·{\mathbf{X}}_{k,l}^{\left(1\right)}+{\widehat{\mathbf{t}}}_0,k=1,\dots ,n_{id,u}l=1,\dots ,n_{id,v}.\hfill \end{array}$$

(17)

The Euclidean distance

$$\begin{array}{cc}\hfill d_{k,l}& =\left|\right|{\mathbf{X}}_{tr,k,l}^{\left(1\right)}-{\mathbf{X}}_{k,l}^{\left(2\right)}\left|\right|\hfill \end{array}$$

(18)

between each transformed point ${\mathbf{X}}_{tr,k,l}^{\left(1\right)}$ belonging to epoch $E_1$ and its corresponding point ${\mathbf{X}}_{k,l}^{\left(2\right)}$ belonging to epoch $E_2$ provides the basis for the decision of whether a surface point is included in the consensus set of model-conforming points. This distance is compared to the error tolerance as follows:

$$\begin{array}{cc}\hfill {ϵ}_{k,l}& =\tau ·{\sigma }_{d_{k,l}}\hfill \end{array}$$

(19)

which is determined by the accuracy ${\sigma }_{d_{k,l}}$ of the computed distance $d_{k,l}$, arising from variance covariance propagation as well as a factor $\tau$, which needs to be chosen appropriately (see below). According to Equation (19), the error tolerance is assessed for each point pair individually since the discretised surface points are not available with a homogeneous precision [30]. When

$$\begin{array}{c}\hfill d_{k,l}\le {ϵ}_{k,l}\end{array}$$

(20)

is fulfilled, the corresponding point pair is assigned to the consensus set ${\mathbb{S}}_{CS}$.

This procedure is repeated either until the number $n_{CS}$ of points fulfilling Equation (20) is larger than or equal the minimum size $n_{min}$ of the consensus set or until the maximum number of iterations $i_{max}$ is reached. It is worth noting that it is not possible to distinguish distortions from gross errors in the data during this step. For this reason, it is essential to successfully remove gross errors a priori using standard techniques, for example, when estimating best-fitting B-spline surfaces.

The successful application of the developed strategy requires the definition of three parameters:

• The parameter $\tau$ in Equation (19) has to be chosen in such a way that only model-compliant points are included in ${\mathbb{S}}_{CS}$. When $\tau$ is chosen too small, points will be erroneously identified as outliers, whereas when $\tau$ is chosen too large, outliers may be included in the consensus set and, thus, the estimation result will become biased. The influence of the parameter $\tau$ on the results is investigated in Section 3.2.
• The minimum number $n_{min}$ of model-compliant points required to accept the current ${\mathbb{S}}_{CS}$ is estimated from the expected amount of gross errors $n_{err}$ ($0\le n_{err}\le 1$) contained in the data set. When $n_{err}$ is known, the optimal $n_{min}$ can be determined according to the following [41]:

  $$\begin{array}{cc}\hfill n_{min}& =(1-n_{err})·n_{id}.\hfill \end{array}$$

  (21)
• The maximum number of iterations $i_{max}$ is determined by the desired probability P that a solution without outliers is found [41]:

  $$\begin{array}{cc}\hfill i_{max}& =\frac{log(1-P)}{log(1-{(1-n_{err})}^{n_{sub}})}.\hfill \end{array}$$

  (22)

Finally, after the iteration is completed, the point pairs of the largest consensus set are used to estimate the transformation parameters ${\widehat{\mathbf{t}}}_{CS}$ and ${\widehat{\mathbf{R}}}_{CS}$ in the extended Gauß–Markov model (14) and (15).

With a final statistical global test, comparing the respective a priori variance factor ${\sigma }_0^2$ with the corresponding a posteriori variance factor ${\widehat{\sigma }}_0^2$, the presence of outliers can be excluded with a confidence probability of $1-{\alpha }_G$, with ${\alpha }_G$ being the error probability [42]:

$$\mathrm{N}\mathrm{u}\mathrm{l}\mathrm{l} \mathrm{h}\mathrm{y}\mathrm{p}\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{s} H_{0,G}:E\left\{{\widehat{\sigma }}_0^2\right\}={\sigma }_0^2$$

(23)

$$\mathrm{A}\mathrm{l}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e} \mathrm{h}\mathrm{y}\mathrm{p}\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{s} H_{A,G}:E\left\{{\widehat{\sigma }}_0^2\right\}>{\sigma }_0^2.$$

(24)

The corresponding test variable $T_G$ is F-distributed with the degrees of freedom f and ∞:

$$T_G=\frac{{\widehat{\sigma }}_0^2}{{\sigma }_0^2}\sim F_{f,\infty }.$$

(25)

If $T_G$ lies outside the associated quantile and, thus, the null hypothesis (23) has to be rejected, a further RANSAC iteration should be performed [30]. (It is worth noting that in this formulation of the global test the presence of errors in the stochastic model is excluded).

2.1.6. Statistically Based Localisation of Distortions

Having successfully conducted the steps described in the previous section, estimated parameters of the rigid body movement as well as a consensus set ${\mathbb{S}}_{CS}$ are available. The latter is a set of $n_{CS}$ point pairs that is free from outliers w.r.t to the estimated rigid body movement with the confidence probability of the global test $1-{\alpha }_G$. Therefore, these points can already be assigned to the non-distorted region. However, depending on the parameter $\tau$ to be chosen, the consensus set does usually not cover the entire non-distorted region of the B-spline surface. Hence, no exact statement can be made about the position and the extent of the distorted regions. Therefore, the consensus set is stepwise and statistically ensured extended by individual point pairs. With the extended consensus set, the parameters of the rigid body movement can be estimated with a higher redundancy and, thus, more precisely.

The starting point of the localisation is the set ${\mathbb{S}}_T$ of points to be tested, which includes all data points that are not assigned to the consensus set at the beginning. Using a forward strategy (see below for information about the order in which the points are considered), these points are then sequentially added to the group of points that is already detected as non-distorted. Based on this extended set, the parameters of the rigid body movement are re-estimated. Point pairs that are detected as outliers within this re-estimation are assigned to the distorted region of the surface, whereas points that do not significantly change the estimated transformation parameters are allocated to the set of non-distorted points, the extended consensus set ${\mathbb{S}}_{CS,ex}$.

In order to implement the outlier detection, the classical Gauß–Markov model is extended by an additional parameter vector $\widehat{\nabla }$ that allows for the estimation of possible outliers [42]. This extension enables the simultaneous testing of several observations with regard to gross errors, a strategy that is desirable in the context of estimating rigid body movements, as a testing of single coordinates is not meaningful [30].

The linearised Gauß–Markov model to estimate the parameters of the rigid body movement (cf. Equations (14) and (15)) is the null hypothesis of the statistical test performed [42]:

$$\begin{array}{cc}\hfill H_0& :E\left\{\mathbf{l}\right\}=\mathbf{l}+\mathbf{e}=\mathbf{A}{\widehat{\vartheta }}_R.\hfill \end{array}$$

(26)

Compared to Equations (14) and (15), the corresponding data points ${\mathbf{X}}^{\left(1\right)}$ and ${\mathbf{X}}^{\left(2\right)}$ are summarised in the observation vector $\mathbf{l}$, and the transformation parameters $\mathbf{R}$ and $\mathbf{t}$ to be estimated are combined in the vector of unknowns ${\vartheta }_R$ in Equation (26). The design matrix $\mathbf{A}$ describes the functional relationship between observations and unknowns.

Using the additional parameter vector $\widehat{\nabla }$, containing the outliers to be estimated, and the corresponding design matrix $\mathbf{B}$, being a sparse matrix with numbers of ones at the entries of the suspected observations, the alternative hypothesis can be formulated as follows:

$$\begin{array}{cc}\hfill H_A& :E\left\{\mathbf{l}\right\}=\mathbf{l}+{\mathbf{e}}^{\prime }=\mathbf{A}{\widehat{\vartheta }}_R^{\prime }+\mathbf{B}\widehat{\nabla }.\hfill \end{array}$$

(27)

Due to the model extension, the parameter vector ${\widehat{\vartheta }}_R$ of the rigid body movement as well as the residual vector $\mathbf{e}$ are modified (indicated by the apostrophe). The dimensions of $\widehat{\nabla }$ and $\mathbf{B}$ are determined by the number $n_a$ of suspected observations. The null hypothesis (26) corresponds to the equivalent null hypothesis $H_0:E\{\nabla \}=0$. Hence, it is examined whether the estimated gross errors differ significantly from zero. The test variable

$$T=\frac{{\widehat{\nabla }}^T{\mathbf{Q}}_{\widehat{\nabla }\widehat{\nabla }}^{-1}\widehat{\nabla }}{n_a·{\widehat{\sigma }}_0^{2^{\prime }}} \sim F_{n_a,f-n_a},$$

(28)

used to investigate the validity of the null hypothesis, follows the F-distribution $F_{n_a,f-n_a}$, the degrees of freedom of which are determined by the number $n_a$ of suspected observations as well as the redundancy f of the initial estimation problem [42]. The computation of T requires the determination of the following quantities:

$$\begin{array}{cc}\hfill {\mathbf{Q}}_{\epsilon \epsilon }& ={\mathbf{Q}}_{\mathbf{ll}}^{-1}-\mathbf{A}{\left({\mathbf{A}}^T{\mathbf{Q}}_{\mathbf{ll}}^{-1}\mathbf{A}\right)}^{-1}\mathbf{A}\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill -\epsilon & ={\mathbf{Q}}_{\epsilon \epsilon }{\mathbf{Q}}_{\mathbf{ll}}^{-1}\mathbf{l}\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill {\mathbf{Q}}_{\widehat{\nabla }\widehat{\nabla }}& ={\left({\mathbf{B}}^T{\mathbf{Q}}_{\mathbf{ll}}^{-1}{\mathbf{Q}}_{\epsilon \epsilon }{\mathbf{Q}}_{\mathbf{ll}}^{-1}\mathbf{B}\right)}^{-1}\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill \widehat{\nabla }& =-{\mathbf{Q}}_{\widehat{\nabla }\widehat{\nabla }}{\mathbf{B}}^T{\mathbf{Q}}_{\mathbf{ll}}^{-1}\epsilon \hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill {\mathsf{\Omega }}^{\prime }& ={\epsilon }^T{\mathbf{Q}}_{\mathbf{ll}}^{-1}\epsilon -{\widehat{\nabla }}^T{\mathbf{Q}}_{\widehat{\nabla }\widehat{\nabla }}^{-1}\widehat{\nabla }\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill {\widehat{\sigma }}_0^{2^{\prime }}& =\frac{{\mathsf{\Omega }}^{\prime }}{f-n_a}.\hfill \end{array}$$

(34)

In the equations above, ${\mathbf{Q}}_{\epsilon \epsilon }$ is the cofactor matrix of the vector of residuals $\epsilon$, ${\mathbf{Q}}_{\widehat{\nabla }\widehat{\nabla }}$ the cofactor matrix of the additional parameter vector $\widehat{\nabla }$ and ${\mathsf{\Omega }}^{\prime }$ the corrected sum of squared residuals.

With the strategy described above, the stability of single point pairs ${\mathbf{X}}_{k,l}^{(1,2)}$ with respect to distortions can be evaluated in a statistically ensured way (single point test).

In order to obtain an areal character of the localisation, alternatively, a predefined domain of the surface can be included in the outlier detection. A straightforward way to define this domain is to use the parametric neighbourhood $[u_{k-a},u_{k+a}]\times [v_{l-a},v_{l+a}]$ of the investigated point pair ${\mathbf{X}}_{k,l}^{(1,2)}={\widehat{\mathbf{S}}}^{(1,2)}(u_k,v_l)$ (see Figure 5 for an example with $a=1$). The parameter a to be chosen determines the size of the neighbourhood. Depending on the choice of a,

$$\begin{array}{cc}\hfill n_a& =3{(2a+1)}^2\hfill \end{array}$$

(35)

observations are examined for gross errors [30]. Naturally, Equation (35) does only hold when ${\mathbf{X}}_{k,l}^{(1,2)}$ does not lie at the surface’s boundary. In that case, $n_a$ reduces accordingly.

Having determined the neighbourhood of the point pair to be tested, all point pairs lying in the domain $[u_{k-a},u_{k+a}]\times [v_{l-a},v_{l+a}]$ are temporarily included in the extended consensus set ${\mathbb{S}}_{CS,ex}$ by adjusting the matrix $\mathbf{B}$ (cf. Equation (24)) accordingly. (It is worth noting that previous test decisions on neighbouring points are not taken into account during this step: neighbouring points that have already been allocated either to the distorted or to the non-distorted regions as well as points that have not been under investigation yet are equally treated and, thus, are involved in the formation of the neighbourhood.) If the comparison of the resulting test variable $T_{k,l}$ (cf. Equation (28)) with the corresponding quantile of the Fisher distribution supports the null hypothesis, the existence of gross errors can be ruled out with the respective confidence probability $1-\alpha$. As a consequence, the point pair under investigation is added to the extended consensus set and is considered non-distorted in the subsequent outlier tests. Otherwise, if the null hypothesis has to be rejected, a distortion of the surface in the domain that is under investigation must be assumed. In this case, the extended consensus set is not modified. With the removal of the point pair ${\mathbf{X}}_{k,l}^{(1,2)}$ from ${\mathbb{S}}_T$, an iteration step of the procedure is completed and the next point pair is investigated. The procedure is continued until all point pairs are removed from ${\mathbb{S}}_T$. At the end, the extended consensus set is used to estimate the final parameters of the rigid body movement. Additionally, a completing global test can be used to check the consistency of the detected non-distorted point pairs.

The order in which the point pairs are investigated by means of the strategy described above, is determined by their degree of consistency with the determined rigid body movement. For this purpose, all points ${\mathbf{X}}_{k,l}^{\left(1\right)}$ to be tested are transformed using the transformation parameters ${\widehat{\mathbf{R}}}_{CS,ex}$ and ${\widehat{\mathbf{t}}}_{CS,ex}$, determined by means of the (extended) consensus set. The transformed points ${\mathbf{X}}_{tr,k,l}^{\left(1\right)}$ are then compared with their correspondences ${\mathbf{X}}_{k,l}^{\left(2\right)}$ by computing the Euclidean distances $d_{k,l}$ (cf. Equation (18)). Since a small $d_{k,l}$ indicates the consensus of the corresponding point pair with the estimated transformation model, of all the point pairs to be tested, the one with the smallest $d_{k,l}$ is examined first [30]. The strategy to localise the distortions is summarised in the schematic sketch in Figure 6.

2.1.7. Regularisation of the System of Equations

By taking into account the accuracies of the sampled surface points, the system of equations for estimating the transformation parameters becomes singular as soon as more surface points are sampled than control points have been estimated (cf. Section 2.1.3). However, in order to achieve a (quasi)-continuous statement regarding the position and extent of the distortions, the number of sampled surface points should considerably exceed the number of control points [30]. The associated singularities of the system of equations must be taken into account accordingly. In [30], three possibilities to deal with the singularities—neglecting the correlations between the estimated surface points, use of the pseudoinverse and regularisation of the VCM’s main diagonal—are investigated and compared. As the second strategy outperforms the first and the third one w.r.t the correctness of the results achieved, the pseudoinverse is used in this contribution if not stated otherwise.

2.2. Data Sets under Investigation

2.2.1. Data Simulation

The B-spline-based strategy to estimate rigid body movements and to simultaneously detect non-distorted regions introduced in Section 2.1 is applied to a variety of simulated data sets in this study, the simulation process of which is demonstrated by means of an example data set. Due to the use of simulated data, the obtained results can be compared with nominal surfaces and nominal transformation parameters and, thus, the developed strategy can be directly validated.

The starting point of the data simulation is the cubic B-spline surface with $(n+1)·(m+1)=7·9=63$ control points presented in Figure 7. This surface is considered the reference surface in the remainder of this section.

The scanning process of this object during the first measuring epoch $E_1$ is emulated by regularly sampling the B-spline surface along the parameter lines, resulting in 10,000 surface points with a spatial resolution of 4–5 mm. Afterwards, the surface points are superimposed by normally distributed measuring noise $\mathbf{n}\sim \mathcal{N}(\mathbf{0},{\mathsf{\Sigma }}_{\mathbf{nn}})$, with ${\mathsf{\Sigma }}_{\mathbf{nn}}$ being the VCM of the measuring noise. In these initial studies, the measuring noise is modelled to be uncorrelated, with a standard deviation of ${\sigma }_{n_x}={\sigma }_{n_y}={\sigma }_{n_z}=\frac{1}{3}$ mm. The choice of these specific values characterising the data points’ precision leads to an average point error of $1$ mm [30].

In order to simulate the point cloud of the second epoch $E_2$, the reference surface is distorted by moving one or more control points. This approach exploits the locality of B-spline surfaces (cf. Section 2.1.1), meaning that the movement of a single control point changes only a local part of the surface. For example, moving control point ${\mathbf{P}}_{3,3}$ (red triangle in Figure 7) by 1.2 cm upwards results in a maximum distortion of the B-spline surface of up to $6$ mm. Figure 8 presents the distances between corresponding points on the reference surface and on the distorted surface in the parameter space (cf. Equation (1) for the relationship between Cartesian coordinates and surface parameters). As can be seen, only parts of the parameter space are influenced by the control point’s movement (the non-distorted parts are not coloured). In most of the distorted area, the distances between corresponding surface points are very small (<1 mm, coloured in green). Compared to the simulated measurement noise, these deviations w.r.t to the reference surface are not significant. The largest deformation occurs in the middle of the distorted part (coloured in red), with a maximum distortion of up to 6 mm.

Afterwards, the distorted B-spline surface is subjected to a similarity transform with neglected scale factor, using the B-splines’ property of invariance w.r.t a similarity transform. For all generated data sets, the following transformations and rotations are chosen:

$${\mathbf{t}}^T=\left[\begin{array}{ccc}0.300\mathrm{m}& 0.600\mathrm{m}& 0.000\mathrm{m}\end{array}\right]$$

(36)

$${\left[\begin{array}{ccc}\omega & \varphi & \kappa \end{array}\right]}^T={\left[\begin{array}{ccc}35.000\mathrm{g}\mathrm{o}\mathrm{n}& 0.000\mathrm{g}\mathrm{o}\mathrm{n}& -10.000\mathrm{g}\mathrm{o}\mathrm{n}\end{array}\right]}^T$$

(37)

The control points of the distorted and transformed B-spline surface, therefore, result from the following:

$$\begin{array}{cc}\hfill {\mathbf{P}}_{ij}^{\left(2\right)}& =\mathbf{R}·{\mathbf{P}}_{ij}^{\left(1\right)*}+\mathbf{t},\hfill \end{array}$$

(38)

with the asterisk indicating the local distortion w.r.t the reference surface.

Finally, a noisy point cloud consisting of 10,000 data points is created by sampling the distorted and transformed B-spline surface and by subsequently adding white noise as described above.

The results of the data simulation are two noisy point clouds describing an object in two measuring epochs that is subject to rigid body movements and superimposed local distortions.

2.2.2. Introduction of the Simulated Data Sets

Figure 9 represents the distortions of the five data sets investigated in this study. The denotation of the data sets contains the percentage of distortion (proportion of the distorted area related to the entire area). In addition, Appendix A

Table A1. Simulated data sets, the maximum distortions $d_{max}$ of which are generated by a shift $d_{\mathbf{P}}$ of one ore more control points. The denotation of the data sets contains the percentage of distortion ${\epsilon }_d$ [30].

Name	Number of Shifted ${\mathbf{P}}_{\mathit{ij}}$	${\mathit{d}}_{\mathbf{P}}$ [mm]	${\mathit{\epsilon }}_{\mathit{d}}[\%]$	${\mathit{d}}_{\mathit{max}}$ [mm]
V0	0	0.0	0	0.0
V15	2	6.0	15	3.9
V18	5	12.0	18	8.1
V20	2	12.0	21	7.8
V40	7	12.0	39	8.1
V65	12	12.0	66	8.1

lists the data sets as well as some of their properties.

As can be seen, the data sets generated differ in the extent of the distorted area (controlled by the number of shifted control points), in the magnitude of the distortion and in the location of the distortions. In order to generate data sets V15 (Figure 9, top left) and V20 (Figure 9, top right), two control points lying in the middle of the control net are shifted. Due to the larger amount of the shift, the distortion in data set V20 has a larger extension and a larger magnitude than the distortion in data set V15. Data set V18 is generated by shifting five control points that are located at the boundary of the control net. As a consequence, the distortion is also located at the boundary of the parameter space (cf. Figure 9, middle left). Finally, data set V40 (Figure 9, middle right) is generated by combining the distortions of data set V20 and V18, and data set V65 (Figure 9, bottom) is generated by adding further distorted regions to V40.

In addition to the distorted data sets, a non-distorted data set V0 is created in order to investigate how susceptible the developed localisation procedure is w.r.t the incorrect detection of distortions in non-distorted regions.

When investigating the simulated distortions (cf. Figure 9), it is worth noting that not only two, but three regions of the surfaces have to be distinguished: In data sets V15, V20 and V18, there are regions that—due to the locality of the B-spline surfaces—are completely non-distorted (not coloured in Figure 9). However, there are also regions in which the surface is distorted, but the distortions of which are so small (≈0$\mathrm{m}\mathrm{m}$) that the distortions are not significant w.r.t to the measurement noise (coloured in green). Finally, there is a region in which the surface is clearly distorted. For a better distinction of this region, a black line is drawn within the green regions of Figure 9, delineating the regions where the deformation is >1 $\mathrm{m}\mathrm{m}$. In data sets V40 and V65, only minimally distorted and distorted regions exist.

3. Results

3.1. Constructing Corresponding Surface Points

The estimation of approximating B-spline surfaces is a well-studied field of research, especially when only the positions of the control points are estimated, as is the case in this study. Hence, the estimation of B-spline surface is not the focus in this contribution. Having determined approximating B-spline surfaces, they are discretised as described in Section 2.1.3. In order to examine the influence of possible singularities arising from the discretisation, different levels of discretisation are investigated in the following sections (cf.

Table 1. Overview of the discretisations [30].

Discretisation	${\mathit{n}}_{\mathit{u}}·{\mathit{n}}_{\mathit{v}}$	${\mathit{n}}_{\mathit{id}}$	Property of ${\mathbf{\Sigma }}_{\widehat{\mathit{l}}\widehat{\mathit{l}}}$
I-V0	$7·9$	63	regular
I-V15	$7·9$	63	regular
I-V20	$7·9$	63	regular
I-V18	$7·9$	63	regular
I-V40	$7·9$	63	regular
I-V65	$7·9$	63	regular
II-V40	$14·18$	252	singular
II-V65	$14·18$	252	singular
III-V40	$28·36$	567	singular
III-V65	$28·36$	567	singular

). With $(n+1)·(m+1)=7·9=63$ estimated control points, all point sets yielded by a discretisation of $n_u·n_v=7·9$ result in regular estimation problems.

3.2. Investigation of the Free Parameters’ Influence on the RANSAC-Based Estimation of Rigid Body Movements

As pointed out, the results of the RANSAC method described in Section 2.1.5 strongly depend on the free parameters to be chosen, especially the parameters $\tau$ (cf. Equation (19)), a (cf. Equation (35)) and $n_{min}$ (cf. Equation (21)).

In order to investigate the influence of these parameters,

Table 2. Influence of selected RANSAC parameters on the number of iterations $\#i$ and the size of the consensus set $n_{CS}$. The results are averaged over 10 applications. ’-’: no solution is found. Top: regular system of equations, bottom: singular system of equations [30].

Discretisation	$\mathit{\tau }$	a	${\mathit{n}}_{\mathbf{min}}$	$\overline{\#\mathit{i}}$	$\overline{{\mathit{n}}_{\mathbf{CS}}}$
I-V0	1	0	12	57	13
	2	0	12	2	27
I-V40	1	0	12	-	-
	2	0	12	20	15
	2	1	16	40	19
	3	1	16	11	23
I-V65	2	0	12	549	12
	3	1	16	214	17
II-V40	2	0	12	37	31
	2	1	40	316	46
	3	0	16	220	37
	3	1	30	684	44
II-V65	2	0	12	86	20
	2	1	20	810	23
	3	1	16	1495	28
III-V40	2	0	12	310	18
	2	1	16	1015	32
	3	0	12	-	-
III-V65	2	0	12	892	16
	3	0	12	-	-

shows the results of the RANSAC strategy (number of iterations $\#i$ and number of point pairs in the final consensus set $n_{CS}$) under variation of the critical parameters for selected discretisations. The results shown are averaged over 10 applications of the strategy and are separated according to the regularity of the resulting system of equations (top: regular system of equations, bottom: singular system of equations) [30].

The choice of $\tau =1$ is investigated for the undistorted data set I-V0 as well as for data set I-V40. Both applications reveal that this choice is too pessimistic: Although a very small consensus set can be found after relatively many iterations for I-V0, no solution is found even after hundreds of iterations for I-V40. When using $\tau =2$ in combination with a regular system of equations, the error threshold (19) is increased and, as a consequence, the number of required iterations is clearly reduced both for data set I-V0 and I-V40. The choice of $\tau =2$ also results in a sufficiently large consensus set for data set I-V65, but requires a very large number of iterations. The further increase to $\tau =3$ leads to a further reduction of the iterations and a larger consensus set. Values of $\tau >3$ are not considered at this point: taking into account the 3-$\sigma$-rule (cf. [1]), it must be expected that these choices are not realistic.

An opposite behaviour is apparent when considering the data sets causing a singular system of equations: The increase from $\tau =2$ to $\tau =3$ leads to a significant increase of the required iterations for almost all examined data sets. Although the pure RANSAC procedure converges after a decreased number of iterations, the resulting consensus sets do not pass the final global test, making further iteration steps necessary. This effect is intensified by a reduced level of discretisation as well as by an increasing amount of distortions. As a consequence, the use of $\tau =3$ does not at all result in a consensus set for which the final global test is successful when considering data sets III-V40 and III-V65. Furthermore, it is noticeable that—although the overall number of points increases significantly by increasing the discretisation density—the size of the consensus set does not.

When a and, thus, the size of the neighbourhood considered in the outlier test are increased, the minimum size $n_{min}$ of the consensus set changes [30]. This behaviour is due to the fact that the maximum number of observations that can be included in the outlier test is influenced by the redundancy of the initial model [42]. Because of the increase of $n_{min}$, the number of required iterations and the size of the resulting consensus set also increase.

The investigations reveal that a threshold determination with $\tau =2$ and $\tau =3$ delivers a sufficiently large consensus set for the subsequent deformation analysis when considering a discretisation with a $7·9$ grid. With $\tau =3$ both a reduction of the required iterations and an enlargement of the final consensus set can be achieved [30]. However, falsifications of the consensus set due to point pairs lying in the slightly deformed area are to be expected more frequently for $\tau =3$. As the subsequent global test rules out the existence of these gross errors with a confidence probability of $1-{\alpha }_G$, nevertheless, $\tau =3$ will be used in the following investigations when considering a regular system of equations. However, the use of $\tau =3$ is not expedient when considering an increased discretisation density. Hence, $\tau =2$ is used in these cases.

Finally, it is worth noting that not only the parameters of the RANSAC algorithm influence the results of the deformation analysis, but also the choice of the significance levels 1-$\alpha$ and 1-${\alpha }_G$: the larger the significance level, the smaller the probability of type I errors, whereas the probability of type II errors can be reduced by appropriately decreasing the significance level. Despite this relationship, the significance levels are not considered tuning constants in this contribution as these parameters are usually fixed or chosen w.r.t the hazard that is monitored.

3.3. Deformation Analysis with a Regular System of Equations

Six discretisations listed in Table 1 (I-V0, ⋯, I-V65) result in a regular system of equations. These data sets are used in this section to demonstrate the general applicability of the developed strategy.

3.3.1. Deformation Analysis for Data Set V0 (No Distortions)

Applying the RANSAC-based estimation of the rigid body movements to the undistorted data set V0 results in the estimated transformation parameters listed in

Table 3. Results of the estimated rigid body movement (data set V0) [30].

${\mathit{t}}_{\mathit{x}}$ [m]	0.300	${\mathit{\sigma }}_{{\mathit{t}}_{\mathit{x}}}$ [mm]	0.023
$t_y$ [m]	0.600	${\sigma }_{t_y}$ [mm]	0.015
$t_z$ [m]	0.000	${\sigma }_{t_z}$ [mm]	0.026
$\omega$ [gon]	35.002	${\sigma }_{\omega }$ [mgon]	3.816
$\varphi$ [gon]	0.008	${\sigma }_{\varphi }$ [mgon]	4.392
$\kappa$ [gon]	−9.997	${\sigma }_{\kappa }$ [mgon]	4.279

. Obviously, the rigid body movement can be estimated very accurately and very precisely for data sets without an additional distortion.

When applying the localisation strategy to data set V0, the majority of point pairs is correctly determined to be non-distorted, but also a few type I errors occur.

Table 4. Results of the localisation for selected point pairs (data set V0). Type I errors are printed in bold [30].

	a	$\mathit{d}\left[\mathrm{mm}\right]$	${\mathit{n}}_{\mathbf{CS}}$	T	${\mathit{F}}_{99}$
${\mathbf{X}}_{0,1}^{(1,2)}$	0	0.206	26	1.17	4.08
	1	0.206	20	1.63	2.36
${\mathbf{X}}_{2,5}^{(1,2)}$	0	0.180	54	4.82	3.91
	1	0.180	47	1.23	1.90
${\mathbf{X}}_{3,4}^{(1,2)}$	0	0.111	26	1.08	4.08
	1	0.111	21	2.23	2.13
${\mathbf{X}}_{6,8}^{(1,2)}$	0	0.717	62	3.48	3.89
	1	0.717	56	1.62	2.30

lists the testing results for selected point pairs of data set V0 in dependence of the parameter a, specifying the neighbourhood used during the testing. In addition to the test variable T and the corresponding 99%-quantile $F_{99}$, the size $n_{CS}$ of the consensus set and the distance d between the respective point pair is given. Type I errors are printed bold. As shown by the example ${\mathbf{X}}_{2,5}^{(1,2)}$, a type I error can be eliminated by considering the direct neighbourhood instead of the single point during the test. However, the exact opposite happens in the case of point ${\mathbf{X}}_{3,4}^{(1,2)}$: being detected as non-distorted in the single point test, it is incorrectly determined as distorted when the direct neighbourhood is taken into account. Thus, the occurrence of type I errors cannot be eliminated by defining a suitable neighbourhood. In any case, tests with a = 1 lead to a homogenisation of the calculated test variables [30]. Remarkably, type I errors occur exclusively in the interior of the surface, whereas edge (e.g., ${\mathbf{X}}_{0,1}^{(1,2)}$) and corner points (e.g., ${\mathbf{X}}_{6,8}^{(1,2)}$) are correctly detected to be non-distorted [30].

3.3.2. Deformation Analysis of Data Sets V15 and V20 (Distortions in the Middle of the Surface)

Table 5. Results of the estimated rigid body movement after having extended ${\mathbb{S}}_{CS}$ (data set V15, $a=0$, $n_{CS}=25$, $n_{CS,ex}=45$) [30].

${\mathit{t}}_{\mathit{x}}$ [m]	0.300	${\mathit{\sigma }}_{{\mathit{t}}_{\mathit{x}}}$ [mm]	0.026
$t_y$ [m]	0.600	${\sigma }_{t_y}$ [mm]	0.017
$t_z$ [m]	0.000	${\sigma }_{t_z}$ [mm]	0.028
$\omega$ [gon]	35.000	${\sigma }_{\omega }$ [mgon]	4.211
$\varphi$ [gon]	0.010	${\sigma }_{\varphi }$ [mgon]	4.912
$\kappa$ [gon]	−9.997	${\sigma }_{\kappa }$ [mgon]	4.848

summarises the final results of the estimated rigid body movements for data set V15 when using the single point test ($a=0$). The initial consensus set of $n_{CS}$ is extended during the localisation procedure by 20 point pairs ($n_{CS,ex}=45$). The resulting estimated parameters of the rigid body movement differ only minimally from those of the non-distorted data set V0 (cf. Table 3). Noticeable differences are only visible in the standard deviations of the estimated parameters, which are slightly larger for data set V15 than for data set V0. The effect is amplified when the direct neighbourhood is included in the tests ($a=1$,

Table 6. Results of the estimated rigid body movement after having extended ${\mathbb{S}}_{CS}$ (data set V15, $a=1$, $n_{CS,val}=14$, $n_{CS,ex}=28$) [30].

${\mathit{t}}_{\mathit{x}}$ [m]	0.300	${\mathit{\sigma }}_{{\mathit{t}}_{\mathit{x}}}$ [mm]	0.037
$t_y$ [m]	0.600	${\sigma }_{t_y}$ [mm]	0.027
$t_z$ [m]	0.000	${\sigma }_{t_z}$ [mm]	0.045
$\omega$ [gon]	34.994	${\sigma }_{\omega }$ [mgon]	6.880
$\varphi$ [gon]	0.000	${\sigma }_{\varphi }$ [mgon]	9.491
$\kappa$ [gon]	−10.008	${\sigma }_{\kappa }$ [mgon]	7.951

). As the initial consensus set is validated by means of individual tests with $a=1$ before performing the localisation (cf. Section 2.1.6), the validated consensus set is smaller than the initial consensus set for $a=0$ ($n_{CS,val}=14$). A similar behaviour can be seen in the sizes of the extended consensus sets: When single point tests are performed, ${\mathbb{S}}_{CS,ex}$ is clearly larger ($n_{CS,ex}=45$)—and, thus, the distorted region significantly smaller—than when the direct neighbourhood is considered during the testing ($n_{CS,ex}=28$). This effect is also evident from Figure 10, presenting for all point pairs the computed test values T of the localisation procedure and the resulting test decision both for $a=0$ (left) and for $a=1$ (right). However, both scenarios ($a=0$ and $a=1$) have in common that, regardless of the choice of the initial ${\mathbb{S}}_{CS}$, the null hypothesis (23) of the global test is accepted after having extended the consensus set.

The estimated transformation parameters achieved by means of data set V20 (analogously to V15 summarised in the Appendix B

Table A2. Results of the estimated rigid body movement after having extended ${\mathbb{S}}_{CS}$ (data set V20, $a=0$, $n_{CS}=21$, $n_{CS,ex}=48$) [30].

${\mathit{t}}_{\mathit{x}}$ [m]	0.300	${\mathit{\sigma }}_{{\mathit{t}}_{\mathit{x}}}$ [mm]	0.031
$t_y$ [m]	0.600	${\sigma }_{t_y}$ [mm]	0.021
$t_z$ [m]	0.000	${\sigma }_{t_z}$ [mm]	0.034
$\omega$ [gon]	35.002	${\sigma }_{\omega }$ [mgon]	5.022
$\varphi$ [gon]	0.020	${\sigma }_{\varphi }$ [mgon]	5.873
$\kappa$ [gon]	−9.999	${\sigma }_{\kappa }$ [mgon]	5.887

and

Table A3. Results of the estimated rigid body movement after having extended ${\mathbb{S}}_{CS}$ (data set V20, $a=1$, $n_{CS,val}=11$, $n_{CS,ex}=20$) [30].

${\mathit{t}}_{\mathit{x}}$ [m]	0.300	${\mathit{\sigma }}_{{\mathit{t}}_{\mathit{x}}}$ [mm]	0.043
$t_y$ [m]	0.600	${\sigma }_{t_y}$ [mm]	0.035
$t_z$ [m]	0.000	${\sigma }_{t_z}$ [mm]	0.072
$\omega$ [gon]	35.002	${\sigma }_{\omega }$ [mgon]	8.789
$\varphi$ [gon]	0.018	${\sigma }_{\varphi }$ [mgon]	17.087
$\kappa$ [gon]	−10.008	${\sigma }_{\kappa }$ [mgon]	10.616

) are almost identical to those yielded by data set V15. The distortions’ dimension and/or magnitude, however, seem to have an influence on the standard deviations of the estimated parameters, as they are larger for data set V20 than for data set V15. Moreover, when using single point tests for data set V20 (Table A2), an unexpected behaviour can be observed: although every null hypothesis of the localisation procedure (26) for point pairs added to the extended consensus set is accepted, the null hypothesis (23) of the final global test has to be rejected. When changing the initial consensus set, it can further be observed that this behaviour does not always occur, but seems to strongly depend on the choice of the initial consensus set. A detailed investigation reveals that point pairs in the transition area between distorted and non-distorted regions are the reason for this behaviour: when those point pairs are assigned to the initial consensus set, the null hypothesis of the global test has to be rejected, whereas it is accepted when only point pairs clearly lying in the non-distorted region are contained in the initial consensus set. Furthermore, it can be observed that increasing the size of the initial consensus set also leads to an increased success of the final global test.

3.3.3. Deformation Analysis of Data Set V18 (Distortions at the Boundary of the Surface)

The results of the deformation analysis of data set V18 fully support the results of V15 and V20: the transformation parameters can be estimated very accurately and very precisely (cf. Appendix B

Table A4. Results of the estimated rigid body movement after having extended ${\mathbb{S}}_{CS}$ (data set V18, $a=0$, $n_{CS}=23$, $n_{CS,ex}=47$) [30].

${\mathit{t}}_{\mathit{x}}$ [m]	0.300	${\mathit{\sigma }}_{{\mathit{t}}_{\mathit{x}}}$ [mm]	0.031
$t_y$ [m]	0.600	${\sigma }_{t_y}$ [mm]	0.019
$t_z$ [m]	0.000	${\sigma }_{t_z}$ [mm]	0.038
$\omega$ [gon]	34.995	${\sigma }_{\omega }$ [mgon]	5.206
$\varphi$ [gon]	0.011	${\sigma }_{\varphi }$ [mgon]	5.986
$\kappa$ [gon]	−10.000	${\sigma }_{\kappa }$ [mgon]	4.973

and

Table A5. Results of the estimated rigid body movement after having extended ${\mathbb{S}}_{CS}$ (data set V18, $a=1$, $n_{CS,val}=17$, $n_{CS,ex}=34$) [30].

${\mathit{t}}_{\mathit{x}}$ [m]	0.300	${\mathit{\sigma }}_{{\mathit{t}}_{\mathit{x}}}$ [mm]	0.038
$t_y$ [m]	0.600	${\sigma }_{t_y}$ [mm]	0.024
$t_z$ [m]	0.000	${\sigma }_{t_z}$ [mm]	0.056
$\omega$ [gon]	34.997	${\sigma }_{\omega }$ [mgon]	7.309
$\varphi$ [gon]	0.094	${\sigma }_{\varphi }$ [mgon]	7.898
$\kappa$ [gon]	−10.003	${\sigma }_{\kappa }$ [mgon]	5.586

), and the distortions can also be successfully detected. When using the direct neighbourhood in the localisation procedure, the region that is detected to be distorted is considerably larger than when single point tests are performed (cf. Figure 11).

3.3.4. Deformation Analysis of Data Sets V40 and V65 (Distortions in the Middle and at the Boundary of the Surface)

As for the data sets considered above, the estimated parameters of the rigid body movement (cf.

Table 7. Results of the estimated rigid body movement after having extended ${\mathbb{S}}_{CS}$ (data set V40, $a=0$, $n_{CS}=19$, $n_{CS,ex}=28$) [30].

${\mathit{t}}_{\mathit{x}}$ [m]	0.300	${\mathit{\sigma }}_{{\mathit{t}}_{\mathit{x}}}$ [mm]	0.041
$t_y$ [m]	0.600	${\sigma }_{t_y}$ [mm]	0.028
$t_z$ [m]	0.000	${\sigma }_{t_z}$ [mm]	0.045
$\omega$ [gon]	35.007	${\sigma }_{\omega }$ [mgon]	5.785
$\varphi$ [gon]	-0.004	${\sigma }_{\varphi }$ [mgon]	8.032
$\kappa$ [gon]	−10.025	${\sigma }_{\kappa }$ [mgon]	6.860

for $a=0$ and

Table 8. Results of the estimated rigid body movement after having extended ${\mathbb{S}}_{CS}$ (data set V40, $a=1$, $n_{CS,val}=7$, $n_{CS,ex}=7$) [30].

${\mathit{t}}_{\mathit{x}}$ [m]	0.300	${\mathit{\sigma }}_{{\mathit{t}}_{\mathit{x}}}$ [mm]	0.141
$t_y$ [m]	0.600	${\sigma }_{t_y}$ [mm]	0.097
$t_z$ [m]	0.000	${\sigma }_{t_z}$ [mm]	0.118
$\omega$ [gon]	35.014	${\sigma }_{\omega }$ [mgon]	12.867
$\varphi$ [gon]	-0.014	${\sigma }_{\varphi }$ [mgon]	20.105
$\kappa$ [gon]	−9.989	${\sigma }_{\kappa }$ [mgon]	21.812

for $a=1$) are very close to the nominal values for data set V40.

However, the precision of the estimated parameters decreases further, both for the single point test and when considering the direct neighbourhood during the testing. Again, the region detected as distorted becomes clearly larger when the direct neighbourhood is considered than when a single point test is performed (cf. Appendix B Figure A1). As for data set V20, a dependency of the initial consensus set’s choice on the results can also be seen for data set V40: point pairs lying in the transition area between distorted and non-distorted region of the surface are sometimes assigned to the distorted region and sometimes to the non-distorted region when the points in the initial consensus set are varied. As a consequence, the null hypothesis of the final global test has to be rejected in the latter case. Again, this behaviour can be reduced by increasing the size of the initial consensus set.

Due to the large amount of distortions, only investigations with $a=0$ are possible for data set V65. The further increase in the amount of distortions compared to the previously investigated data sets results again in a rejection of the null hypothesis of the final global test. However, due to the small non-distorted regions, it is not possible to control the success of this data set’s deformation analysis by increasing the size of the initial consensus set.

3.4. Deformation Analysis with a Singular System of Equations

Using data set III-V40 as an example, the results of a deformation analysis with a singular system of equations are shortly presented in this section. Figure 12 (left) shows the results of the localisation procedure for data set III-V40 when performing a single point test. Comparing this result to Figure A1 (left), presenting the equivalent result for I-V40, the large amount of type II errors lying in the transition zone between the two distorted regions is striking. However, as already described above, the number of type II errors can be clearly reduced when the size $n_{CS}$ of the consensus set is increased. This effect is clearly visible in Figure 12 (right), showing the results of the localisation procedure when the consensus set is larger than in Figure 12 (left). Obviously, the separation of the two distorted regions succeeds much more satisfactorily with the larger consensus set.

4. Discussion

The developed strategy to estimate rigid body movements and simultaneously detect distorted regions is applied to a variety of synthetic data sets, the distortions of which vary in dimension and shape. The results presented in Section 3 reveal that—for justified selected free parameters of the algorithm—the estimated parameters of the rigid body movement are very close to the nominal values for all investigated data sets. Naturally, the larger the distorted region, the smaller the precision of the estimated parameters: due to the reduced number of point pairs in the non-distorted region, the redundancy and, consequently, the precision of the estimated parameters decrease.

The localisation of the distorted areas on the basis of the data set V0 reveals that the occurrence of type I errors cannot be completely avoided. However, since type I errors usually do not occur in groups, most of them can be detected by investigating their direct neighbourhoods: if their entire neighbourhood is located in the non-distorted region, it can be assumed that a type I error exists. Furthermore, with a large confidence probability $1-\alpha$, the probability of the occurrence of type I errors can be regulated.

The results of the data sets that indeed represent distorted objects show that clearly distorted regions can be successfully detected in all investigated scenarios, independently of their location on the object. The only difficulty is to precisely distinguish the actual distorted regions from the non-distorted ones. Particularly in the minimally distorted regions—usually, the transition region between distorted and non-distorted regions—, type II errors occur when performing single points tests: Point pairs lying in the minimally distorted region are erroneously allocated to the non-distorted region. As a consequence, the final global test (${\alpha }_G=5\%$) fails. This failure of the final global test can be reduced by considering the direct neighbourhood during the localisation procedure ($a=1$) rather than conducting single point tests ($a=0$). Using this strategy, the minimally distorted region can be more successfully detected.

Hence, the success of the deformation analysis depends on a suitable choice of the free parameters: Satisfactory results are obtained with $a=1$ for scenarios with relatively large non-distorted regions (e.g., V15, V20 and V18), as is obvious when comparing the detected distorted regions with the nominal ones. For data set V18 (results: Figure 11 (right), nominal: Figure 9 (middle left)), the non-distorted region is almost perfectly delimited from the distorted one. Only the discretised point $\widehat{\mathbf{S}}(1.00,0.63)$ actually belongs to the distorted region, but is assigned to the non-distorted region. However, as the actual distortion in this point is almost zero, the impact of this erroneous allocation is minimal. Similar results are achieved for data set V15 (results: Figure 10 (right), nominal: Figure 9 (top right)). Here, even seven points are erroneously allocated to the non-distorted region, all of them lying on the surface’s boundary. However, the points’ positions are not the reason for the wrong allocation, but the minimal magnitude of the distortion.

The investigations using data sets V40 and V65 show that it is possible to satisfactorily estimate rigid body movements, even when there are no non-distorted regions at all. Alternatively, there must be regions in which the distortions are so small that they nevertheless are allocated to the consensus set. This behaviour can be observed particularly well by means of data set V40 (results: Figure A1 (right), nominal: Figure 9 (middle right)): As for data sets V18 and V15, few discretised points that are only minimally distorted are erroneously allocated to the non-distorted region. Nevertheless, the estimated parameters of the rigid body movement are very close to the nominal parameters. However, the larger the distorted regions are, the more emphasis has to be placed on the definition of the initial consensus set. The larger this consensus set is, the more reliably type II errors can be avoided, as the initial estimation of the rigid body movement succeeds more accurately.

Furthermore, the results of data set V40 show that also for data sets with large distorted regions, more accurate localisation results are achieved when considering the direct neighbourhood during the localisation ($a=1$) than when conducting single point tests ($a=0$). However, for data sets in which the distorted regions occupy the majority of the surface (data set V65), a localisation with $a=1$ is no longer possible, and single point tests are the only option. In this case, the behaviour of single point tests that points that lie in the minimally distorted region are increasingly detected as non-distorted must be accepted. To take advantage of the localisation using $a=1$, it has to be ensured during data acquisition that a sufficiently large non-distorted area of the objects is captured.

Singularities arising from a decreased discretisation level can be successfully handled by using the pseudoinverse in all subsequent analysis steps. This approach allows to significantly increase the resolution of the results and, thus, to detect even small-scale distortions.

5. Conclusions

In this contribution, a strategy is introduced that allows for the estimation of rigid body movements of an object and for the simultaneous detection of distorted regions based on laser scanning data. The starting points of the developed strategy are best-fitting B-spline surfaces that approximate the point clouds acquired in at least two measuring epochs. Using the surface parameters u and v, point correspondences on these B-spline surfaces are defined. Based on these point correspondences, rigid body movements between the surfaces are estimated. When the rigid body movements are superimposed by a local distortion, the estimated parameters are falsified. In order to detect these gross errors w.r.t to the rigid body movement, a robust RANSAC-based strategy to estimate rigid body movements is implemented. The result of this step is an initial solution for the rigid body movement as well as the consensus set ${\mathbb{S}}_{CS}$ that is a subset of the surfaces’ non-distorted region. Hence, the consensus set is extended in a second step using statistical hypotheses tests within a forward strategy: Point pairs that do not significantly change the estimated transformation parameters are also allocated to the non-distorted region. With the final non-distorted region, the parameters of the rigid body movement are re-estimated.

Developed as an S2S comparison, the measurement noise and the amount of data is reduced to a large extent during the deformation analysis. The results obtained are approved by means of statistical tests, allowing for the significance assessment of the determined deformations.

The developed strategy is applied to a variety of simulated data sets, the local distortions of which vary in shape, extent and magnitude. For all investigated data sets, the parameters of the rigid body movement can be estimated very accurately and very precisely. The detection of the distorted regions succeeds in all data sets when the local neighbourhood is considered during the localisation procedure. When single point tests are conducted, type II errors occur as points in the transition area are erroneously allocated to the non-distorted region.

Point pairs lying in the distorted region can then be directly used to determine displacement vectors. Alternatively, the parameter grids can be investigated within a strain analysis, resulting in an areal deformation statement rather than in a point-wise one. This strain analysis based on the B-spline parameter grid is the content of future investigations, just like the application of the developed strategy to measured data sets. In contrast to the study presented here, in which the stochastic model of the data is known due to the use of simulated data, the determination of a realistic stochastic model of the data is of superior importance in this future research.