Variances of Cylinder Parameters Fitted to Range Data

Industrial pipelines are frequently scanned with 3D imaging systems (e.g., LADAR) and cylinders are fitted to the collected data. Then, the fitted as-built model is compared with the as-designed model. Meaningful comparison between the two models requires estimates of uncertainties of fitted model parameters. In this paper, the formulas for variances of cylinder parameters fitted with Nonlinear Least Squares to a point cloud acquired from one scanning position are derived. Two different error functions used in minimization are discussed: the orthogonal and the directional function. Derived formulas explain how some uncertainty components are propagated from measured ranges to fitted cylinder parameters.


Introduction
Digital representation of an object's surface is important in many engineering applications. Examples include: manufacturing, modular construction, and preservation of historical artifacts. Frequently, the data needed for such representation are collected with 3D imaging systems. The systems are non-contact, lineof-sight instruments which provide range images ( , ) r ϑ ϕ of the object's surface, where r denotes the distance (measured as a time of flight) from an instrument to a point on a surface, and ϑ and ϕ are the elevation and the azimuth angles to that point. The scanning process is very fast and current systems can collect datasets containing hundreds of thousands of points within a few seconds [1]. Points collected from a surface may be used to build a 3D geometrical model of a scanned object. For example, a cylinder is frequently modeled in projects where industrial pipelines are scanned [2][3][4][5]. Frequently, Nonlinear Least Square (NLS) methods [6][7][8][9][10] are used to fit a model to its corresponding point cloud. An important issue is how to propagate an instrument uncertainty (for example: uncertainty of measured ranges r) to the uncertainties of fitted cylinder parameters. These uncertainties are essential in a quality control process when the as-built model (derived from the point cloud) is checked against the as-designed (theoretical) model [11]. The tolerances used in the design process and the uncertainties of the fitted model parameters are both needed to accept the as-built model. For example, the uncertainties of the fitted coefficients of the cylinder centerline are needed to check, with a given confidence, if two pipes are parallel or to check the plumbness of a cylinder column.
It is commonly accepted by manufacturers and users of 3D imaging systems that the uncertainty of a j-th measured point ( , , ) j X Y Z P is caused mainly by an uncertainty of the measured range r j ; an impact of angular uncertainty can be neglected and both the elevation and the azimuth angles of P j are treated as noise-free parameters. Under this assumption, the closed formulas for variances (i.e., squared standard deviations) of cylinder parameters fitted to data acquired from one scanning location are provided. These variances are an important component in evaluating measurement uncertainty. The choice of an error function for NLS minimization is a delicate matter. It was shown that the orthogonal error function, which is commonly used in commercial software, yields results that are in disagreement with experiment when a plane is fitted to a point cloud acquired in some experimental conditions [12]. Another error function, the directional error function, gives results that are, on average, in a much better agreement with experiment. Therefore, in this paper, closed formulas for variances of fitted cylinder parameters are derived for both orthogonal and directional error functions. Both error functions depend on four fitting variables: three angles { } , , ϑ ϕ ρ which determine a cylinder centerline orientation in space, and the distance L of a centerline to the origin of a coordinate system.

Variances of Fitted Cylinder Parameters
on a surface of an infinitely long regular cylinder of radius R satisfy the following equation where the cylinder centerline M is ( , ) ϑ ϕ w is a unit vector in the direction of the axis, parameterized by two angles: the elevation ϑ and the azimuth ϕ , and × stands for a cross product of two vectors. Therefore, the Cartesian coordinates of ( , ) L is any fixed point on the axis and q is a number defining a location of arbitrary point on the line M. It is convenient to choose L in such a way that its length L = L is the shortest [5]. This happens when L and w are perpendicular and their dot product is zero, 0 ⋅ = L w . In this setting, L is the distance of the cylinder centerline M from the origin of the coordinate system. The origin and the orientation of the coordinate system (in which experimental points P j are acquired) are defined by the location and the orientation of a scanner. Vector L may be parameterized as where a unit vector l lies on a plane perpendicular to w and The unit vectors u and v, together with the vector w, form the right-handed orthonormal basis {u, v, w} The angle ρ in Eq. (4b) defines two coordinates of the unit vector l on the 2D plane perpendicular to w, as shown in Fig. 1 where r j is a range and a unit vector p j is equal to ( , ) [cos cos , cos sin ,sin ] It is also convenient to introduce a vector d j Now, using the above introduced vector d j together with Eqs. (4) and (5), Eq. (1) can be expressed as Thus, in this parameterization, a regular cylinder of radius R is defined by four parameters: three angles , , ϑ ϕ ρ and the distance L (in this and the next two sections, a radius R is treated as a known constant; later in Section V this restriction is removed for the orthogonal fitting). Due to imperfections in range measurement, experimental points P j do not lie exactly on a cylinder surface. Therefore, a model of a cylinder has to be fitted to the experimental dataset P where E j is the squared distance between the experimental point P j and its corresponding theoretical point on a cylinder surface. Different definitions of the theoretical point yield different error functions. In this paper we study two error functions: the orthogonal error function E O and the directional error function E D , as explained in Figs. 2 and 3 and in the next two sections. The elements of a gradient of the error function E ∇ can be calculated as The knowledge of closed formulas for a gradient E ∇ allows application of efficient minimization algorithms to find the fitted parameters { } * * * * , , ,L ϑ ϕ ρ . Their values depend solely on the experimental dataset P {N} , i.e., The variances of the fitted cylinder parameters * var( ) ϑ , * var( ) ϕ , * var( ) ρ , and * var( ) L may be calculated following the same general approach developed originally for fitting a sphere to range data [13]. Here, similarly as in the previous study, only range uncertainty is considered and the bearings ( ) , j j ϑ ϕ of acquired j-th point P j are treated as noise-free parameters. For other types of instruments, for example coordinate measuring machines (CMM), the above assumption may not be valid and the formulas for variances of fitted cylinder parameters developed in this paper may not be applicable. In addition, it is assumed that measured range r j is not correlated with r k for any j ≠ k. When both assumptions are valid, the variances of the fitted cylinder parameters may be estimated by applying to Eq. (13) the uncertainty propagation formula [14] ( ) For every j, the sensitivity vector    (15) can be evaluated by solving the following 4x4 system of linear equations where the vector V j is defined by while H is 4x4 Hessian matrix which elements are given by

Orthogonal Fitting
When fitting a cylinder using the orthogonal error function (see Fig. 2) the theoretical point corresponding to the experimental point P j is defined as the orthogonal projection of P j on a cylinder surface. Thus, Eq. (11) takes the form where Q j is a distance of point P j to the cylinder axis M, which can be evaluated from Eq. (10) ( , , , , ) ( sin Defining the scalar functions a j and b j by means of Eqs.

Directional Fitting
When fitting a cylinder using the directional error function, see Fig. 3, the theoretical point T j corresponding to the experimental point P j is defined as an intersection of a ray originating from the instrument and passing through P j with the cylinder surface where the function ( ) , , , , j j t L ϑ ϕ ρ P has its value close to one when optimization variables are close to their best-fit values and the theoretical points T j satisfy Eq. (1) of the cylinder ( ) A ray need not always intersect a cylinder surface -a similar situation occurs in the directional fitting of a sphere [15]. Whether there is an actual intersection or not can be checked by calculating the distance h j between the cylinder centerline M given by Eq.
(2) and a ray passing through P j ; if h j < R, then there is an intersection and the theoretical point on a cylinder surface is defined by Eq. (23). Otherwise, if h j ≥ R, there is no intersection and the theoretical point T j is defined as a point on a cylinder surface which has the smallest distance to a ray passing through P j . Assuming that a ray is not parallel to the cylinder centerline The j-th term in Eq. (11) corresponding to the directional error function E D is a squared distance between the experimental point P j and the theoretical point T j , see Fig. 3 ( 1) ( , , , , ) From Eq. (28), the j-th term in Eq. (12) for the components of the gradient of the directional error function

Discussion
So far, only four fitting variables { } , , ,L ϑ ϕ ρ were used and the radius of a cylinder R was treated as a known constant. When the radius is unknown and needs to be fitted, the orthogonal error function E O in Eq. (19) may be still effectively minimized in four dimensional space { } , , ,L ϑ ϕ ρ [5]. This is possible because the radius R enters Eq. (19) in a quadratic term. Since a gradient of the error function must be zero when the minimum of E O is reached, the radius R can be evaluated from the condition Substituting for E O Eq. (19), the radius R may be calculated as where Q j is calculated in Eq. (22). Equation (31) has a simple geometric interpretation: the fitted radius * R is the algebraic mean of distances of all experimental points P j from the fitted cylinder centerline defined by { }  (14) using the formulas derived in Section III for the orthogonal fitting. Unfortunately, the above strategy for fitting the unknown cylinder radius R cannot be applied to the directional fitting because a j-th term in the directional error function E D depends on R in a more complex way and a convenient separation of fitting variables as in Eq. (31) is not possible.
Closed formulas for variances derived in this paper are applicable to a cylinder fitted to a point cloud acquired from only one scanning location. When two or more datasets acquired from different locations are registered together, the formulas presented here cannot be used. This limitation has its root in the way range uncertainties are propagated into uncertainties of fitted parameters: the registration process makes relatively straightforward formulas much more complicated, with variances explicitly depending on the registration parameters. It remains an open question whether fitting a model to one large point cloud consisting of many registered datasets has an advantage over fitting a model to many unregistered clouds.
The difference between the orthogonal and the directional fitting becomes more important for point clouds containing many points acquired at large angle of incidence (i.e., AOI -the angle between a ray coming through the data point P j and a normal to the cylinder surface at the theoretical point T j defined in Eq. (23)). Scanning industrial pipelines (where a length of a scanned pipe may be larger than the fitting variable L) may yield a dataset in which AOIs are varying in a large range. Then, the orthogonal fitting may generate a cylinder centerline which is biased and incorrect. In order to quantify the expected differences between both error functions, further research involving computer simulations and field measurements is needed.

Appendix
This Appendix contains expressions needed for evaluating the directional error function E D , defined in Eq. (28). The individual j-th term of E D depends on the distance h j between the cylinder centerline M and a ray passing through data point P j . Let us first consider the case when h j , given by Eq. (26), is less than the cylinder radius R. Eqs. (27) and (24) lead to the following equation for t j where t j depends on P j and parameters determining M, i.e., , , ϑ ϕ ρ , and L. By squaring both sides of Eq. (A1) and using the definitions of a j and b j from Eqs. (21a) and (21b), the following quadratic equation for t j is obtained This quadratic equation has two solutions. From a definition of the theoretical point T j in Eq. (23), it follows that only t j > 0 has a physical meaning. Two different experimental settings for collecting data points P {N} need to be distinguished: a) when a cylinder (for example, a pipe) is scanned from the outside; and b) when scanning is done from inside (for example, in a tunnel or inside a containment vessel for a nuclear reactor).
In the first case, L > R and b j defined in Eq. (21b) is always positive (a j from Eq. (21a) is always positive, regardless of a scanner location). Then, from the two real solutions of Eq. (A2), a physical http://dx.doi.org/10.6028/jres.117.015 meaning exists only for the smaller one which corresponds to the intersection of a ray with a cylinder surface closer to a scanner location (the origin of coordinates), as shown in Fig. 3. It is useful to define an auxiliary function ( , , , , ) Once t j is known, a j-th term in the directional error function E D can be calculated for the case when h j < R according to Eq. (28). When h j ≥ R , scanning can be done only from outside and a ray passing through P j does not intersect with a cylinder surface; then quadratic Eq. (A2) does not have a real solution and t j cannot be determined. In such a case, a point G j on a ray is sought which has the smallest distance to the cylinder centerline M, that is The actual value of s j can be evaluated by finding a minimum of a squared distance of G j to M, where w is given by Eq.  The distance h j defined in Eq. (26) together with s j provided in Eq. (A10) are used to calculate a j-th term in the directional error function E D when h j ≥ R, according to Eq. (28).