Measurement uncertainty assessment for virtual assembly

Virtual assembly (VA) is a method for datum definition and quality prediction of assemblies considering local form deviations of relevant geometries. Point clouds of measured objects are registered in order to recreate the objects’ hypothetical physical assembly state. By VA, the geometrical verification becomes more accurate and, thus, increasingly function oriented. The VA algorithm is a nonlinear, constrained derivate of the Gaussian best fit algorithm, where outlier points strongly influence the registration result. In order to assess the robustness of the developed algorithm, the propagation of measurement uncertainties through the nonlinear transformation due to VA is studied. The work compares selected propagation methods distinguished from their levels of abstraction. The results reveal larger propagated uncertainties by VA compared to the unconstrained Gaussian best fit. 1 Current trends in dimensional metrology and state-of-the-art datum definition and uncertainty assessment As quality demands on products increase, tolerance specifications for parts become more and more complex. With these challenging geometrical specifications, verification algorithms are required that represent the geometrical system more precisely. According to Nielsen (2003), in the last few decades, dimensional tolerances shrank due to improved manufacturing systems. However, the form deviations could not be reduced by the same extent. Therefore, their consideration should be intensified. A main deficit in the current International Organization for Standardization (ISO) standard for datum definition, ISO 5459 (Deutsches Institut für Normung e.V., 2011), is the lack of consideration of local form deviations for datum features. A datum feature is defined as a “real (non-ideal) integral feature used for establishing a single datum” (Deutsches Institut für Normung e.V., 2017, p. 2). Datum systems composed of three datum features mathematically define a coordinate system. This allows the definition of tolerance zones for extrinsic tolerances (Weißgerber and Keller, 2014). About 80 % of all measurement tasks require datum systems, so a further function-oriented datum system definition has a strong impact on geometrical verification. Hence, an assessment of the uncertainty for datum systems is of broad interest. Figure 1 shows a datum definition, where three perpendicular associated planes are considered in a nested approach. The primary datum constrains 3 degrees of freedom (DOF), the secondary datum 2 DOF and the tertiary datum 1 DOF (Gröger, 2015). 1.1 Concept of the virtual assembly In this paper, measurement data of physical objects are gathered from measurements using industrial computed tomography (CT). Registration is the action of aligning a data set relatively to another according to a datum definition in a common coordinate system. Virtual assembly (VA) comprises the consideration of local form deviations in the datum system computation. As shown in Fig. 1a, through VA, the physical workpiece contact is simulated by computing the contact points. The registration for VA is mathematically stated as an optimization problem, as introduced in Weißgerber and Keller (2014). In the following, matrices are marked as boldface capital, vectors in boldface italic, and scalar Published by Copernicus Publications on behalf of the AMA Association for Sensor Technology. 102 M. Kaufmann et al.: Measurement uncertainty assessment for virtual assembly Figure 1. Datum definition by nested registration, using associated planes (a), registration approach according to the default case in the current standard, (b) and registration approach according to virtual assembly (c). values in roman formatting. The signed distances dsig,i of i ∈ 1. . .N , i ∈ N, corresponding pairs of points { p1,i,p2,i } , with p1,i ∈ P 1 and p2,i ∈ P 2 determine the clearance between the surfaces to register. P 1 and P 2 are point sets of surfaces 1 and 2, respectively, as presented in Fig. 1b and c. The objective function f of the optimization problem is as follows: f ( tx, ty, tz,φ,θ,ψ )


Current trends in dimensional metrology and state-of-the-art datum definition and uncertainty assessment
As quality demands on products increase, tolerance specifications for parts become more and more complex. With these challenging geometrical specifications, verification algorithms are required that represent the geometrical system more precisely. According to Nielsen (2003), in the last few decades, dimensional tolerances shrank due to improved manufacturing systems. However, the form deviations could not be reduced by the same extent. Therefore, their consideration should be intensified. A main deficit in the current International Organization for Standardization (ISO) standard for datum definition, ISO 5459 (Deutsches Institut für Normung e.V., 2011), is the lack of consideration of local form deviations for datum features. A datum feature is defined as a "real (non-ideal) integral feature used for establishing a single datum" (Deutsches Institut für Normung e.V., 2017, p. 2). Datum systems composed of three datum features mathematically define a coordinate system. This allows the definition of tolerance zones for extrinsic tolerances (Weißgerber and Keller, 2014). About 80 % of all measurement tasks re-quire datum systems, so a further function-oriented datum system definition has a strong impact on geometrical verification. Hence, an assessment of the uncertainty for datum systems is of broad interest. Figure 1 shows a datum definition, where three perpendicular associated planes are considered in a nested approach. The primary datum constrains 3 degrees of freedom (DOF), the secondary datum 2 DOF and the tertiary datum 1 DOF (Gröger, 2015).

Concept of the virtual assembly
In this paper, measurement data of physical objects are gathered from measurements using industrial computed tomography (CT). Registration is the action of aligning a data set relatively to another according to a datum definition in a common coordinate system. Virtual assembly (VA) comprises the consideration of local form deviations in the datum system computation. As shown in Fig. 1a, through VA, the physical workpiece contact is simulated by computing the contact points. The registration for VA is mathematically stated as an optimization problem, as introduced in Weißgerber and Keller (2014). In the following, matrices are marked as boldface capital, vectors in boldface italic, and scalar values in roman formatting. The signed distances d sig,i of i ∈ 1. . .N , i ∈ N + , corresponding pairs of points p 1,i , p 2,i , with p 1,i ∈ P 1 and p 2,i ∈ P 2 determine the clearance between the surfaces to register. P 1 and P 2 are point sets of surfaces 1 and 2, respectively, as presented in Fig. 1b and c.
The objective function f of the optimization problem is as follows: where the sum of the squared signed distances d sig is minimized. The optimization variables t x , t y , and t z determine the translation, and φ, θ, and ψ are the Euler angles of the rigid transformation of the point set to register P 2 to the fixed point set P 1 . The avoidance of a physical intersection of parts is formulated as the omission of surface intersection, where P 1 ∩ P 2 = ∅. This optimization constraint can be either formulated as a hard constraint, as implemented in this work, allowing d sig,i ≥ 0 only, or by a soft constraint by introducing a penalty term in order to permit small intersections. However, the constrained optimization is formulated as an unconstrained, differentiable optimization problem by using the Lagrange multiplier method (Griva et al., 2009). Thus, a Lagrange function L, according to the following: with λ as the Lagrange multiplier and g (M) as sum of squared negative signed distances d sig,neg,j is introduced, where j ∈ 1; N neg , j ∈ N + , and N neg = D sig,neg is the cardinality of the set of negative signed distances D sig,neg = d sig | d sig < 0 for a particular transformation M t x , t y , t z , φ, θ, ψ . The first-order optimality condition can be stated as ∇L = 0. The transformation of each point p 2,i of the point set P 2 from its initial scene coordinate system into the transformed world coordinate system, defined by the datum system and, in the following, denoted by superscript Tr, is defined as follows: where P Tr 2 is the transformed point set, R is the rotation matrix gathered by matrix multiplication of the individual rotation matrices for φ, θ , and in the named order, T is the translation matrix with its components t x , t y , and t z , and M is the composed transformation matrix. The signed distances d sig,i are determined according to the following: where n i is the normal vector in P 1,i . The optimization problem is nonlinear due to trigonometric functions arising at transformation, as shown in the explicit depiction in matrix notation in Eq. (5) for a point p 2,i p 2,i,x , p 2,i,y , p 2,i,z , as per the following: cos θ · cos ψ − cos θ · sin ψ sin θ tx sin φ · sin θ · cos ψ + cos φ · sin ψ − sin φ · sin θ · sin ψ + cos φ · cos ψ − sin φ · cos θ ty − cos φ · sin θ · cos ψ + sin φ · sin ψ cos φ · sin θ · sin ψ + sin φ · cos ψ cos φ · cos θ tz (5)

Introduction to measurement uncertainty assessment
A complete statement of a measurement result includes the measurement uncertainty. The measurement uncertainty is a nonnegative quantity expressing doubt about the measured value, defined as a "parameter, associated with the result of a measurement, that characterizes the dispersion of the values that could reasonably be attributed to the measurand" (ISO/IEC, 2008b). In ISO/IEC (2008a), the main stages of uncertainty assessment are described as formulation, propagation, and summarizing. During formulation, generally the measurand Y is defined, input quantities X are determined, and the measurement model is established as follows: A main step is the uncertainty estimation for the input quantities. The Guide to the expression of uncertainty in measurement (GUM) (ISO/IEC, 2008b) reveals two methods for uncertainty estimation. With a type A evaluation, statistical quantities from independent, uncorrelated observations are assessed. With a type B evaluation, the uncertainty is estimated based on a priori information, such as specifications, previous measurements, calibration data, or user experience. Since no a priori knowledge is available in this work, the type A evaluation is conducted as described in Sect. 3.1. Uncertainties are separated in systematic and stochastic contributors. Since no calibrated values are available, the systematic error (bias) from a calibrated true value is omitted in this work. During the propagation step, uncertainties assigned to Y are propagated through the measurement model. An overview on propagation methods is given in Sect. 3.2. At the summarizing step, the propagated uncertainty is consolidated into a coverage interval, extended by coverage factor k (ISO/IEC, 2008a).

Aim and scope of this paper
The aim of this work is to analyze the uncertainty propagation for the VA algorithm. Since only few contact points may influence the hard constraint of the optimization problem, a lower robustness compared to existing methods is initially assumed. At the moment, uncertainty propagation is commonly not considered for fitting algorithms such as the VA. However, information gathered from the uncertainty propagation could, on the one hand, be used to claim the robustness of a registration result and, hence, the derived measurements. On the other hand, this allows a reduction in the uncertainty contribution of registration algorithms, for example, by avoiding less certain registration results. As a use case, a linear guide assembly, consisting of a slider mounted to a rail, is assessed, which is shown in Fig. 3. Measurement data were captured using the CT system Werth TomoScope HV 500, with an acceleration voltage of 180 kV, a tube current of 240 mA, an integration time of 500 ms, 1000 projection images, and a resulting voxel size of 0.2 mm.

Description of methods for uncertainty assessment and propagation
The propagation model corresponds to the model function P Tr 2 = M · P 2 outlined in Eq. (3). Each quantity has an associated uncertainty. The initial uncertainties before propagation comprise the transformation uncertainty u T , associated to transformation M, and the point uncertainty U Pt , associated to the point set P 2 . The propagated uncertainty U Tr is associated to the transformed measurement point set P Tr 2 . In the following, multiple modalities for the declaration of the propagated uncertainty are introduced. The matrix U Tr = u Tr x , u Tr y , u Tr z comprises the uncertainty components of the x, y, and z coordinates. The composition of the uncertainty components u Tr x , u Tr y and u Tr z to the combined standard uncertainty u Tr c is defined as follows: which allows the representation of multivariate measurands in a scalar value. For the mathematical formulation of the propagation, the uncertainties are stated in form of their corresponding covariance matrices. All uncertainties arise from normal distributions due to a Gaussian distribution of random variables in the measurement process. In Fig. 2, the graphical decomposition of the propagation model is shown. In the following Sect. 3.1, the assessment of the initial uncertainties, U Pt and u T , is described. In Sect. 3.2, the propagation to determine U Tr and u Tr c is described.

Assessment of the initial uncertainties
The measurement point uncertainty U Pt associated to P 2 is assessed experimentally by 20 repeated CT scans at constant measurement settings. It is determined according to the single point uncertainty approach proposed by Fleßner et al. (2016). For measurement point i ∈ 1. . .N , i ∈ N + , the uncertainty matrix is defined according to the following: The column components u Pt,i,x = u Pt,i,y = u Pt,i,z = u Pt,i are equal, since the used uncertainty assessment approach does not allow the consideration of a probing direction. The uncertainty u Pt,i equals the standard deviation s Pt,i of the normally distributed distances d v,i that are gathered from the 20 scan repetitions v ∈ 1. . .V , with d i as the average distance of all repetitions for each particular point index i, according to the following: In order to calculate the Euclidean norm, as follows: which is referring to the average point p i , corresponding sets of points p i , p v,m are gathered by a k-nearest-neighbors spatial search that identifies the indices m in p v corresponding to each i in P i . P i denotes the point set of average points p i ∈ P i that are the mean coordinates of all repetitions v determined for each i. The uncertainty of transformation is as follows: and it includes the uncertainties of the optimization variables that are stated in Eq. (1). If u T is small, a high repeatability of  registration results can be assumed. The uncertainty u T is assessed by a Monte Carlo propagation, according to ISO/IEC (2008a), as the standard deviations of the optimization variables from K = 10 000 repetitions, where the point set P 2 is slightly transformed with the delta transformation as follows: δT = δt x , δt y , δt z , δφ, δθ, δψ .
The small transformations in δT are random realizations from a population N (0, σ T ) for the translations t x , t y , t z and N (0, σ R ) for rotations φ, θ , and ψ. In this work, based on scientific judgment (type B evaluation according to ISO/IEC, 2008a), the values σ T = 0.1 mm and σ R = 0.001 rad are applied. The mean transformation of all K observations is considered as the true transformation, as follows:

Uncertainty propagation method
For the mathematical formulation of the propagation, the variance-covariance matrix (hereafter just covariance matrix) is introduced. Generally speaking, this quadratic matrix comprises the uncertainties u l,m in form of their variances s 2 l,m = u 2 l,m , where the diagonal entries for l = m represent the variances, and correlations are expressed by the covariances for l = m. The covariances are zero because the input variables are assumed to be uncorrelated, according to Galovska et al. (2012). Eigenvectors and eigenvalues of describe a rotational uncertainty ellipsoid, as shown in Fig. 4. For the generic measurement model Y = f (X), the covariance matrix x is propagated into y by matrix multiplica-tion, according to the following: where J is the Jacobian matrix developed in x, x is the estimate of X, and J is the transposed Jacobian matrix. The derivatives occurring due to the Jacobian matrix J are analogous to a Taylor series expansion. To cover nonlinear propagation functions, higher-order terms of the Taylor series need to be considered. The work of Galovska et al. (2018) gives a further overview on uncertainty propagation methods in the context of virtual measurements of gaps for the automotive body in white. In our work, both the covariance matrices Pt,i for a certain point i, and T gathered from the transformation uncertainty u T , are propagated through the transformation function f (refer to Eq. 3) into Tr (Chong and Mori, 2001;Galovska et al., 2012;Rüschendorf, 2016;Zeier et al., 2012).
Pt,i , and T are formulated as diagonal matrices according to the following: Pt,i , u 2 Pt,i , u 2 Pt,i and (15) The transformed uncertainties u Tr x , u Tr y , and u Tr z correspond to the square root of the diagonal entries in matrix Tr for column and row indices l = m equal to one, two, and three for the x, y, and z components, respectively.
The linear propagation (LP) is valid for small rotations due to the small angle approximation, where sin x ≈ x and cos x ≈ 1 for small values of x. Due to a preceding rough Table 1. Comparison of propagation methods concerning relative difference and normalized computation time in arbitrary units.

Propagation method Linear propagation (LP) Second order propagation (SO) Monte Carlo propagation (MC)
Mean relative difference d rel (standard deviation of d rel ) 2.57 × 10 −3 (1.25 × 10 −3 ) 3 × 10 −6 (7 × 10 −6 ) 0 Normalized computation time referring to LP 1 1.6 14.3 Figure 4. Uncertainty components u Tr x , u Tr y , and u Tr z of the x, y, and z components, respectively (from left to right), and uncertainty ellipsoid for U Pt (small red sphere) and U Tr (large brown ellipsoid) around an arbitrarily selected point. alignment, rotations smaller than 0.01 rad (≈ 0.6 • ) are determined by VA, resulting in relative errors of 1.7 × 10 −5 for the sine function and 5 × 10 −5 for the cosine function. To assess the effect of possible approximation errors, the results for the LP are compared to the propagation results for the second order (SO) and the Monte Carlo (MC) propagation. For this purpose, the propagation methods implemented in the UncLib library in MATLAB R2018b supplied by Wollensack (2020) were used. Figure 3 shows the coordinate system (CS) orientation and the datum of A and B, respectively, that define the assembly. The CS position is centered in the barycenter of the slider. Each point p 2,i in P 2 is rotated around the CS origin p 0 (0 | 0 | 0), with a radius r i , which is defined as the Euclidean distance of the point and origin, according to the following:

Discussion
Due to the contribution of the transformation uncertainty u T to the transformation, the propagated uncertainty U Tr increases with an increasing radius r, which can be related by analyzing the visualizations in Fig. 4, where the components u Tr x , u Tr y , and u Tr z increase towards the corners of the analyzed slider. The uncertainties of rotational transformation parameters cause parabolic and circular patterns for the propagated uncertainties. Here, a large uncertainty u of the rotation in causes a large uncertainty u Tr x . The components u Tr y and u Tr z are affected to a smaller extent by u φ and u ψ , respectively. The uncertainties u T,x , u T,y , and u T,z cause an offset to u Tr x , u Tr y , and u Tr z , respectively. The datum A (see Fig. 3) constraining the x translation shows the largest propagated uncertainty u Tr x of about 0.1 mm. The uncertainty u Tr y of about 0.06 mm, associated to datum B and constraining the y direction, is nearly half as large. The uncertainty u Tr z in the z direction is about 0.01 mm, which matches the initial uncertainty before propagation, since t z is constrained to zero. The fact that u Tr x is about 2 times the magnitude of u Tr y might be an effect of imbalanced point sets because datum B contains about half the points of datum A. Due to a high uncertainty u Tr x , the projection of the uncertainty ellipsoid shown in Fig. 4 in the x-y plane is distorted in the x direction. If no rough pre-alignment is performed before VA, an increase in the uncertainty u Tr c of about a factor of 10 is observed because the optimization result converges poorly for the equal termination criteria. This emphasizes the need for a sufficient pre-alignment. By comparing the uncertainty u Tr c for the VA algorithm to an unconstrained Gaussian best fit registration, an increase in u Tr c of about 50 % to 68 % was observed, as shown in Fig. 6. Hence, the VA algorithm results in approximately twice the measurement uncertainty. The mean relative differences are as follows: For the LP and SO propagation methods referring to the MC propagated uncertainty u Tr c,MC are presented in Table 1. According to ISO/IEC (2008a), the MC method should be  used when a linearization of the propagation function is inadequate and, therefore, is considered as the ground truth method. The propagated distribution is gathered due to propagating a large number of random samples (here 10 6 samples) through the model function. For both the LP and SO propagation, the values of d rel are negligibly small, so that no significant variation in u Tr c was observed. For a median propagated uncertaintyũ Tr of approximately 0.25 mm, the relative difference of 2.57 × 10 −3 for the LP propagation results in an absolute error of approximately ±0.7 µm only. However, relative to the LP propagation, the computational cost for the MC propagation is about 14.3 times higher, while the factor of the SO propagation is merely about 1.6.

Conclusion and outlook
The main contribution to the propagated uncertainty is due to the uncertainty in the transformation parameters u T , which depend on the formulation of the registration problem. Confirming the work of Galovska et al. (2012), it was shown that the number and arrangement of the reference points considered for registration strongly influence the propagated uncertainty. By an unconstrained Gaussian best fit, all points are weighted equally, which reduces the uncertainty compared to the VA algorithm, where certain points are considered in the constraint function. Thus, compared to the unconstrained Gaussian best fit, the uncertainty is nearly doubled. A preceding rough alignment helps to strongly reduce the propagated uncertainty. Due to small transformations remaining for the VA after pre-alignment, the small angle approximation and, hence, a linear propagation model can be sufficiently applied. The propagated uncertainties are relatively large. The main contribution to the propagated uncertainty U Tr is the uncertainty of transformation u T . Large observed values for u T are also assumed to be caused due to large values for the delta transformation δT. Since, in practice, no random shift of the point set to register occurs, for practical measurements the uncertainty of transformation u T and the propagated uncertainty U Tr are assumed to be considerably smaller. As an outlook, the propagation model can be prospectively used to state uncertainties for virtual measurements that are performed on the registered data sets. Here the analysis of the gap flush between the taillight and tailgate for an automotive application by virtual measurements regarding uncertainties is to mention Galovska et al. (2018). Parameter studies on best termination criteria can be performed, aiming to minimize the uncertainty of transformation u T .
The propagation model can be verified by analyzing repeatability studies on virtual measurements of registered assemblies. Moreover, computationally efficient minimum variance estimators, such as Kalman filters can be studied in order to evaluate the preliminary VA registration result during an iteration based on the magnitude of uncertainty. As a further approach, the Jacobian matrix J of the uncertainty propagation can be evaluated during the VA optimization, which is a measure for the sensitivity of the propagated uncertainty towards the transformation variables. By doing so, fewer uncertain registration results can possibly be determined.