Least squares best-fit geometric elements taking into account uncertainty structure

In coordinate metrology, a key activity is fitting a geometric surface to coordinate data. By far the most common fitting criterion has been ordinary least squares (OLS). OLS is appropriate if the covariance matrix associated with the data is a diagonal matrix. However, in recent years much effort has been devoted to developing more realistic uncertainty models for coordinate data and it is now timely to develop fitting algorithms that take into account this richer uncertainty structure. For many measuring systems, the uncertainty characterisation is such that the covariance matrix associated with coordinate data can be represented in a compact, factored form that makes explicit the contribution of random and systematic effects associated with the measuring system. This compact representation enables the least squares fitting algorithm to be implemented using essentially the same computational components as for the OLS fitting scheme.


Introduction
In coordinate metrology, a key activity is fitting a geometric surface to coordinate data.By far the most common fitting criterion has been ordinary least squares (OLS).OLS is appropriate if the covariance matrix associated with the data is a diagonal matrix.However, in recent years much effort has been devoted to developing more realistic uncertainty models for coordinate data and it is now timely to develop fitting algorithms that take into account this richer uncertainty structure.In fact, the effort in developing realistic uncertainty characterisations is only beneficial if the algorithms acting on the coordinate data use this extra information.For many measuring systems, the uncertainty characterisation is such that the covariance matrix associated with coordinate data can be represented in a compact, factored form that makes explicit the contribution of random and systematic effects associated with the measuring system.This compact representation enables the least squares fitting algorithm to be implemented using essentially the same computational components as for the OLS fitting scheme.In addition, the new algorithm can often give insight into the behaviour of the measuring system that gave rise to the measurement data.In particular, the algorithm can be used explicitly to determine uncertainties in coordinate metrology when the coordinate measuring system is used in a comparator mode.
In section 2 we discuss the uncertainty characterisation of measuring systems, while in sections 3 and 4 we discuss, respectively, ordinary least squares (OLS) and generalised least squares (GLS) surface fitting algorithms.Our summary and concluding remarks are given in section 5.
a e-mail: ian.smith@npl.co.uk

Notation
Throughout this paper, we use the following notation to represent the coordinates associated with a set of data points in three dimensions.Given point coordinates

CMM kinematic model
The kinematic error of a coordinate measuring machine (CMM) describes the systematic error behaviour of the CMM arising from non-ideal geometry relating to scale errors, straightness and rotational errors.The error behaviour is defined in terms of the 18 univariate functions d xx , d xy , . .., e zz , six a function of x, six a function of y and six a function of z.These error functions can be combined into a single error function e(x, b) depending on parameters b = (b 1 , . . ., b p ) T .The actual centre of the CMM probe, x * , is related to the measured probe centre x output by the CMM through a model equation of the form where is a random effect.We make the assumption that e is sufficiently smooth so that e(x * , b) can be approximated effectively by e(x, b), i.e., the error function e evaluated at the measured coordinates x.
Estimates of the kinematic error parameters b can be determined experimentally from the measurements of hole or ball plates [9,11], for example.For a ball plate experiment, let y j be the location of the jth ball centre in a fixed frame of reference and let y j,k = T (y j , t k ) = x k,0 + R T (α k )y j DOI: 10.1051/ C Owned by the authors, published by EDP Sciences, 2013 be its location in the kth position of the ball plate, j = 1, . . ., n Y , k = 1, . . ., n T .Here T (y, t) represents a rigid body transformation defined by six parameters t, three translation parameters and three rotation angles.The measurement equations are of the form If the 3 × 1 random effects i are modelled as draws from a Gaussian distribution i ∈ N(0, σ 2 M I), then estimates of the model parameters y j , t k , and a are determined by minimising where a nonlinear least squares problem.If J is the Jacobian matrix of partial derivatives of f i with respect to the optimisation parameters evaluated at the solution, the variance matrix V for the fitted parameters is estimated by V = (J T J) −1 .The variance matrix V b associated with the kinematic error parameters b is a submatrix of V .
In the following, we assume that estimates of the kinematic error parameters, b, have been determined and that, following the implementation of an error correction, b is assigned a multivariate Gaussian distribution.We assume further that e(x, 0) = 0. Similarly, we also assume that the random effects are also described by a Gaussian distribution: Suppose measured data x i , i ∈ I = 1, . . ., m, is generated according to the model defined by where b and i are characterised as in expressions (2).Then, to first order, the 3m × 1 vector of actual probe centres x * I is associated with the multivariate Gaussian distribution where the variance matrix V I is given by 3 Ordinary least squares (OLS) feature assessment Suppose u → s(u, a) defines a parametric curve or surface.The parameters u determine the position of a point on the surface and the parameters a determine the shape and position of the surface.We now assume that set of measured coordinates, x I generated according to the model ( 3) is such that for some u i and some a.
The ordinary least squares (OLS) best-fit surface can be found by minimising where u * i specifies the point s * i = s(u * i , a) on the surface closest to x i and n i is the normal vector at s * i .For standard geometric elements, d(x, a) can be evaluated analytically.For more general surfaces, numerical methods are required [7].Let J be the Jacobian matrix associated with d(x i , a) at the solution a * and let be its QR factorisation [8].Then, to first order, a = where If N is the m × 3m block diagonal matrix storing n T i in the 1 × 3 diagonal blocks, the sensitivity matrix S a of the OLS fitted parameters a with respect to the random effects e I is given by so that the variance matrix associated with the x I is propagated through to the variance matrix associated with a according to (10) The second term on the right is the contribution of the random effects to the variance matrix associated with fitted parameters a and is the term that is usually evaluated in a nonlinear least squares model fit (see, e.g., [5,4]).The first term on the right is the contribution from the parametric errors.The matrix S a S b is the sensitivity of the fitted parameters a j to the parametric errors b k .
To first order the residual error vector is given by d = Q 2 Q T 2 h I with h I as specified in expression (9) and the associated variance matrix is As for V a , if V I is of the form (4), then V d is similarly decomposed into contributions arising from parametric errors and random effects.

Example: effect of scale and squareness on cylindricity assessment
To illustrate this approach, we consider a model of a CMM with scale and squareness errors [2,3]  is associated with N(0, σ 2 I) with σ = 2 × 10 −5 .For this statistical model, the standard uncertainties associated with the point coordinates are of the order of 0.002 mm.A cylinder can be specified by five parameters a = (x 0 , y 0 , α, β, r 0 ) T , where x 0 and y 0 are the coordinates of the point of intersection with the plane z = 0, the rotation angles α and β determine the axis direction and r 0 is the radius.For this calculation, a = (60, 60, 0, 0, 25) T .Table 1 gives the sensitivities of the fitted parameters with respect to the parametric errors.The scale errors have an effect on the location and radius parameters; the xy-squareness parameters b xy has an effect on the y-coordinate of the axis location while the squareness errors b xz and b yz have an effect on the axis direction.Table 2 gives the sensitivities of the residual distances associated with the five points on the middle circle; sensitivities for points on the other two circles are the same.The table shows that only the x-and y-scale errors along with the xysquareness have an effect on the evaluated form error.

Generalised least squares (GLS) feature assessment
The OLS feature assessment corresponds to the maximum likelihood estimate of the feature parameters a for the case V I = σ M I.For more general variance matrices, the maximum likelihood estimate is given by the solution of min a,{ui} where Even for standard geometric features, there is in general no analytical expression for the optimal u * i .However, if V I has the form as in expression (5), where G is a 3m × p matrix made up of 3 × p blocks G i , then the optimal a that solves problem (12) is also the optimal a for the problem min a,b, Thus, the optimal u i defines the distance from the point x i − G i b to the surface u → s(u, a).For standard geometric elements, this distance can be evaluated analytically.Problem (13) is a standard nonlinear least squares problem that can be solved using the Gauss-Newton algorithm in which it is required to evaluate where n i is the 3 × 1 vector of partial derivatives of d(x, a) evaluated at x = x i − G i b.Compared to the OLS fit, the only new information required is the evaluation of n i , a straightforward calculation.If required, the complex-step method can be used to provide an accurate numerical estimate of these derivatives [1,10].For example, if the software that evaluates d(x, a) supports complex arithmetic, then the imaginary part of d(x + ihe, a)/h, where e = (1, 0, 0) T , i = √ −1, and h = 10 −100 , evaluates ∂d/∂x, etc.
The variance matrix Ṽ associated with the fitted parameters according to problem (13) is given by where J and K are defined as in expressions ( 14), evaluated at the solution.The advantage of the GLS formulation ( 12) is that the uncertainty information is used more effectively so that the uncertainties associated with the fitted parameters are smaller that those associated with the OLS fitted parameters as evaluated in expression (10).In terms of J and K, we have Thus V aOLS is the same as V aGLS but with M replaced by I. Since M − I is positive semi-definite, V aOLS − V aGLS must also be positive semi-definite.
Table 3 shows the uncertainties associated with the parameters of a cylinder fitted to data generated according to the model described in section 3.1 using OLS and GLS fits.The uncertainties associated with x 0 and y 0 are significantly smaller for the GLS fit.The eigenvalues of the matrix M −1 are 0.06, 0.12, 1.0, 1.0, 1.0, and 1.0, showing that in a two dimensional subspace, the uncertainty contribution from the systematic scale and squareness effects is much reduced when the GLS fit is used.

Comparator measurements
The advantage of formulation (13), compared to formulation (12), is that the systematic effects b appear explicitly allowing updated estimates of the parameters to be obtained.This advantage opens the way to evaluating the reduction in uncertainties that can be achieved when a CMM is used in comparator mode.
For example, suppose a CMM with scale and squareness errors as described in section 3.1 is used to measure a calibrated ring gauge.We assume that the uncertainty associated with the calibration value of the gauge radius is 0.000 2 mm and the form errors of the gauge are drawn from a Gaussian distribution with standard deviation σ F = 0.000 2 mm [6].This calibration information allows us to derive improved estimates of the x-and y-axis scale errors and the xy-squareness parameters.The uncertainties associated with these estimates are reduced from 2 × 10 −5 to 6.0 × 10 −6 for the x-and y-axis scale errors and 7.3 × 10 −6 for the xy-squareness error.The fourth column in table 4 shows the uncertainties in the fitted parameters for the measurement of a ring gauge using these improved estimates.

Summary and concluding remarks
This paper has been concerned with fitting geometric elements to coordinate data that has an associated uncertainty structure that reflects the systematic effects associated with the measuring system.The uncertainty structure allows us to determine the generalised least squares estimate in a way that is computationally efficient and requires only a small modification to software that implements a standard least squares fit.The approach also allow determines fitted surfaces with associated uncertainties for a CMM used in comparator mode.
and S b is the 3m × p matrix of partial derivatives of e(x i , b) with respect to the error parameters b evaluated at b = 0. Equation (4) decomposes the variance associated with the best estimate x I as the sum of the parametric error and random contributions.If V b has Cholesky decomposition [8] V b = L b L T b , then

b
xx b xy b xz 0 b yy b yz 0 0 b zz   x.The diagonal elements represent the scale errors along the three axis while the non-zero off-diagonal elements are the squareness errors.We assume there are 15 measured points on a cylinder of radius 50 mm, five points on three circles 25 mm apart.We assume that σ M = 0.000 5 mm and that b = (b xx , b yy , b zz , b xy , b xz , b yz ) T

Table 1 .
Sensitivities of cylinder parameters to parametric errors.

Table 2 .
Sensitivities of residual distances to parametric errors.

Table 3 .
Uncertainties associated with fitted parameters.