Variances of Plane Parameters Fitted to Range Data

Formulas for variances of plane parameters fitted with Nonlinear Least Squares to point clouds acquired by 3D imaging systems (e.g., LADAR) are derived. Two different error objective functions used in minimization are discussed: the orthogonal and the directional functions. Comparisons of corresponding formulas suggest the two functions can yield different results when applied to the same dataset.

where • stands for the dot product of two vectors. The Cartesian coordinates of w (ϑ,ϕ) may be written as (2) The absolute value of parameter D is the distance from the plane to the origin of the coordinate system, and D may be expressed as (3) where P 0 is any point on a plane. A plane can be fit to the experimental dataset P {N} = {P j , j = 1,…, N}, where N denotes the number of points; the goal being to calculate the numerical values of the three parameters defining the plane: ϑ, ϕ and D. Within the framework of the Least Squares method, the fit parameters are obtained by minimizing the error function (4) where E j is the distance between the experimental point P j and its corresponding "theoretical point." 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 Fig. 1 and in the next two sections.
It is not surprising that due to nonlinear dependence of the normal vector w on both angles ϑ and ϕ, plane fitting requires nonlinear minimization. However, as is shown in the next two sections, for both error functions E O and E D the distance E j depends linearly on the third parameter D . Therefore, D may be explicitly expressed as a function of both angles (ϑ , ϕ) and P {N} from the condition (5) For any error function E defined by Eq. (4) with E j depending linearly on D, the linear parameter can be expressed as a function of the remaining non-linear parameters and dataset P {N} When the error function E reaches its minimum at [ϑ * , ϕ * , D * ], a gradient of the function has to be zero. This implies that the original 3D search space may be reduced to a 2D space and the error function may be re-written as (7) the minimum of E being located at (8) The variances of the fitted plane parameters var(ϑ * ) and var(ϕ * ) may be calculated following the same general approach developed for fitting a sphere to range data [8], [9]. In the current study, the same assumption is made as in the previous studies: for the 3D imaging systems relevant to this study, the uncertainty in the range measurement is typically much larger than the uncertainty in the angular measurements. Thus, an acquired point P j can be expressed as where r j is a range measured at bearings (ϑ j , ϕ j ) and || P j || = r j . In this approximation the bearings are treated as noise-free control variables and a unit vector p j is defined as Note that 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 plane parameters developed in this paper may not be applicable. In addition, it is assumed that the correlation in the measured ranges r j and r k is negligible for any j ≠ k . When both assumptions are valid, the variances of the fitted plane parameters may be estimated by applying to Eq. (8) the uncertainty propagation formula [2] (10a) (10b) and the covariance may be estimated as (10c) The variance of the third parameter D * may be calculated from the uncertainty propagation formula [2] applied to a general function D (ϑ, ϕ, P {N} ) defined in Eq. (6), where the derivatives of D are calculated at (ϑ * , ϕ * , P {N} ).
The individual sensitivities and (11) may be calculated as in [8] by solving for each j the following 2 × 2 system of linear equations  ( , , ).
The matrix H is the Hessian of the error function These general formulas are now applied to two specific error functions: the orthogonal error function E O and the directional error function E D .

Orthogonal Fitting
For the orthogonal plane fitting (see Fig. 1) the theoretical point O j corresponding to the experimental point P j is defined as the orthogonal projection of P j on a plane. Thus, Eq. (4) takes the form (15) Applying condition (5), Equation (6) can be expressed as (16) where P 0 here is the centroid of all experimental points P {N} . This condition states that the plane fitted with the orthogonal error function has to contain the centroid, P 0 . Defining a scalar product d j of two vectors w and p j given by Eqs. (2)  where r j is a measured range. In this notation, Eq. (7) describing the error function in the reduced 2D search space of angles (ϑ, ϕ) can be written as .
Finally, the elements of vector V j defined in Eq. (13b) can be obtained by differentiating with respect to r j the elements of the gradient ∇E O given by Eqs. (20a) and (20b). Taking into account the definition of vector U j given by Eq. (19b) and the definition of centroid P 0 of all points P j as well the dependence of P j on r j given by Eq. (9a), the elements of vector V j can be evaluated as where (22c) δ j, k is the Kronecker delta, and p k is the unit vector defined in Eq. (9b). First and second order derivates of the vector w (ϑ,ϕ), defined in Eq. (2), are provided in Appendix A, Eqs. (A1-A5). Once the matrix H and vectors V j are known, the sensitivity vectors S j can be calculated for every j by solving the 2 × 2 system of linear Eq. (12). The variances of the fitted angles var(ϑ * ) and var(ϕ * ) and the covariance cov(ϑ * ,ϕ * ) can then be determined from Eqs. (10a, b, c). The variance of the third parameter D * , defined in Eqs. (16-18), can be now evaluated from Eq. (11) using the following equations: where the derivatives of D are calculated at [ϑ * , ϕ * , P {N} ].

Directional Fitting
For the directional plane fitting (see Fig. 1) the theoretical point D j corresponding to the experimental point P j is defined as an intersection of a ray originating from the instrument through P j with the plane where a parameter t j has its value close to 1, and the theoretical points satisfy Eq. (1) of the plane (24b) The distance E j in Eq. (4) , , . the Kronecker delta. Once the matrix H and vectors V j are known, the sensitivity vectors S j can be calculated for every j by solving a 2 × 2 system of linear equations (12). The variances of fitted angles var(ϑ * ) and var(ϕ * ) and the covariance cov(ϑ * ,ϕ * ) can then be determined from Eqs. (10a, b, c). The variance of the third parameter D * , defined in Eq.
Finally, the elements of vector defined in Eq. ( 13b)  Figure 1 shows that the orthogonal plane fitting behaves differently from the directional fitting. The difference between orthogonal and directional fitting depends on the Angle of Incidence (AOI) of the laser beam. For AOI approaching 90º, the optimal value of the orthogonal error function E O (ϑ * , ϕ * , P {N} ) is decreasing, even when the uncertainty of the measured ranges r j is large. The optimal value of the error function is usually interpreted as a gauge of noise level in experimental data (assuming that a right model is fitted to the data). For 3D imaging systems, due to a divergence of a laser beam, range measurements collected for large AOI have large uncertainty. Thus, the behavior of E O is in a sharp contrast with the common experimental observation. The directional fitting is free of this flaw and a residual value of the directional error function E D (ϑ * , ϕ * , P {N} ) correctly estimates a level of noise in the acquired experimental dataset P {N} for any AOI. When E O is minimized, the sensitivities of the large AOI. The flawed sensitivities entered in the Eqs. (10) and (11) will cause an underestimation of the variances of the plane parameters fitted with the orthogonal function for large AOI. For small AOI, the difference between E O and E D is diminishing and both error functions are expected to provide correct estimates for the variances of fitted parameters.

Discussion
Individual sensitivities S j of fitted angles are calculated from Eq. (12). The vector V j on the right hand side of this equation behaves differently for the orthogonal and the directional error function, see Eqs. (D1, 2) and (D7, 8) in the Appendix D. This may cause a much larger variability of S j calculated for the orthogonal fitting and a poorer estimate of variances of fitted parameters than for the directional fitting.
As was already pointed out, a plane fitted by minimizing the orthogonal error function has to contain the centroid P 0 of the acquired points P {N} , see Eq. (16). Directional fitting does not have this constraint. Thus, a minimization of two error functions discussed in this paper may lead to different results.

Conclusions
In this paper we derived formulas for the variances (which are useful for uncertainty analysis) of plane parameters fitted to a dataset acquired with 3D imaging systems. Two error functions were investigated: the orthogonal and the directional error function. Comparison of corresponding formulas suggests the two functions may yield different results when applied to the same range data. However, in order to quantify the anticipated difference, laboratory experiments and computer simulations are needed.

Appendix A
From the definition of w (ϑ , ϕ ) in Eq. (2)

Appendix B
Auxiliary functions defined for calculation of derivatives of D (ϑ , ϕ , P {N} ) defined in Eq. (28) for the directional fitting are given by: Their respective derivatives are: where the derivatives of d j are defined in (A6-A10). , , .