1

 Sensitivities and Variances for Fitted Parameters of Spheres 

 By 

 Christoph Witzgall and Marek Franaszek 

Abstract 

Experimental data have been gathered by applying 3D imaging systems, such as 

LIDAR/LADAR instruments, to spherical objects. This report provides a compilation of the 

statistical and analytical procedures to be used for an evaluation, to be reported separately, of 

two different methods of modeling objects, directional and orthogonal fitting, based on those 

data. Estimating the variances of fitted parameters directly from their sensitivities to data 

perturbations is proposed. Sensitivities are determined by implicit differentiation of error 

gradiens. Detailed descriptions of the directional and orthogonal fitting methods, as applied to 

spheres in a scanning environment, are set forth. In particular, the report furnishes closed-form 

expressions for those derivatives of the respective error functions which are needed for the 

calculation of the parameter sensitivities with respect the full set of control variables. 

Key Words 

3D imaging systems, directional fitting, least squares, LIDAR/LADAR, orthogonal distance 

fitting, parameter variance, parameter sensitivity, sphere, statistical modeling

1 Introduction 

A frequent task is to determine the shape characteristics, size, position and pose of physical 

objects for purposes of identification, location, registration, and calibration of coordinate frames. 

Such tasks are needed, among others, for quality control in manufacturing, determination of “asbuilt” structures, construction automation and site monitoring, e.g., [1, 2]. 

A common approach to these tasks is to acquire 3D coordinates of data points considered to lie 

on the surface of targeted objects. 3D Imaging Systems, which include “line-of-sight” 

LIDAR/LADAR devices, are increasingly used for this purpose. The latter instruments, in 

particular, are capable of fast generation of large amounts of data points or “point clouds”. They 

scan an object by emitting laser pulses and processing return signals in order to determine the 

distance traveled and thus determine the distance or “range” between the instrument and the 

point of impact -- presumably -- on the object. The device keeps track of the “bearings” such as 

azimuth and elevation angle, at which each particular signal was emitted. This process of data 

acquisition suggested the use in this work of polar/spherical (“angle-angle-range”) coordinate 

systems for representing data points. Also, transformation to Cartesian coordinates will introduce 

correlation. 
2

Once the point cloud corresponding to an object has been determined, a computational process is 

required to extract the desired features of the object from this data set. In typical applications, a 

mathematical “model” is specified, based on features characteristic for a class of objects. The 

model is “parameterized”, that is, it is defined with the help of parameters that determine these 

characteristics. Choosing values for these parameters will result in the mathematical description 

of a surface to represent a “virtual object”, which may then be compared to an image of the real 

object as provided by the point cloud. By adjusting the model parameters so that the virtual 

object moves into a location that optimizes the proximity of the object surface to points in the 

point cloud, desired characteristics such as location, pose, size and shape are found within the 

coordinate frame of the point cloud. This permits determining the geometric relationship 

between that object and other objects or features that are also represented in the point cloud 

frame. If this frame registers to an established ground-truth frame, then absolute measurements 

of location, pose and shape can be extracted. 

Such approaches to the modeling of objects of interest within a point cloud may employ the 

powerful “Iterative Closest Point (ICP)” method [3], or the “Hough Transform”, e.g. [4]. Present 

work focuses on the extensively used “Fitting” paradigm, which is based on minimizing a 

specified error function or on maximizing likelihood. The reader may want to consult texts on 

“Statistical Models” such as [5-7]. 

Of particular interest are two least-squares based approaches, “orthogonal” and “directional“ 

fitting. Orthogonal fitting, also referred to as “Orthogonal Distance Regression (ODR)” [8, 9], or 

“Geometric Fitting” [10], is a commonly used and widely commercialized method. In particular, 

publications [10-19] discuss its application to the fitting of spheres or circles. The alternate 

approach, “directional fitting” has been proposed and discussed [20, 21] for data acquired by 

scanning from a single instrument position. Here, the orthogonal (closest Euclidean) distance to 

the virtual object has been replaced by the distance in the direction of the scan by which the data 

point had been acquired. While computational aspects dominate much of this research, our 

interest here is in statistical and metrological issues. 

The thrust of this report is an approach to determining the sensitivities of fitted model 

parameters, in general, and for spherical models, in particular. The report is also preparatory to 

an experimental study of different fitting methods and their statistical evaluation [22]. At issue, 

in particular, is the estimation of derived variances for fitted sphere centers based on specified 

variances for range measurements. In Chapter 2, the general fitting paradigm, based on the 

concept of an error function, is described, along with the general computational formalism for 

calculating parameter sensitivities. These sensitivities will be used to estimate parameter 

variances. In dealing with spherical models, this approach is of necessity more general than the 

common nonlinear least squares approach based on linearization and homoscedacity. A 

comparative discussion of these statistical procedures will be provided in a separate report [23]. 

Chapters 3 and 4 are dedicated, respectively, to orthogonal and directional fitting of spheres in a 

scanning environment. Closed forms of the derivatives, needed for calculating sensitivities, are 

reported. The Appendix will feature detailed derivations of the reported formulas so as to enable 

verification. 
3

2 The Fitting Paradigm 

2.1 Error Function 

Once a parameterized model has been selected, it is natural to ask for parameters that minimize 

the extent to which the point cloud deviates from the resulting virtual object. The hope is that 

such an – at least locally – optimal virtual object provides, within the coordinate frame of the 

data points, an accurate representation of the actual object. Fitting a 3D model of a sphere of a 

Cartesian center C = [X , Y, Z] and known radius R may be accomplished by specifying an 

“error function” 

 

 (2.1.1) ( , , , , , , , , ..., , , , )

1 1 1 2 2 2 n n n E = E X Y Z d φ θ d φ θ d φ θ

where X , denote the model parameters, and the variables Y, Z d i n i i i,φ ,θ , = 1,..., , are 

coordinates of points [ ] i i i i P = d φ θ to be measured. The following discussions, however, 

should not be construed as pertaining only to this special scenario, but rather as representative of 

full generality. In particular, the data may also be Cartesian, or not be coordinates, at all. 

The choice of the error function should be such that it produces only nonnegative values. A 

minimum of zero should indicate a perfect fit. An error function E thus furnishes a model 

description. 

Given an actual data set of measurements 

 [ ], ...,1 , ,

)0( )0( )0( )0( Pi = di φi θi

i = n

the parameter values 

 

)0( )0( )0( X = X , Y = Y , Z = Z

are thus determined by minimizing the expression 

 ( , , , , , , , , ..., , , , )

)0( )0( )0( )0(

2

)0(

2

)0(

2

)0(

1

)0(

1

)0(

1 n n n E = E X Y Z d φ θ d φ θ d φ θ

for the variables X , , given the coordinate values of the data points Y, Z Pi

. 

A common approach to constructing error functions is to assign an individual error 

 ) ( , , , , ,

i i i i i

e = e X Y Z d φ θ

to each data point [ ] i i i i P = d φ θ , and to minimize the sum of squares 
4

(2.1.2) ∑=

=

n

i

i E e

1

2

. 

Both, the orthogonal and directional fitting method, mentioned in the Introduction, are based on 

the Nonlinear Least Squares (NLS) concept [17, 19]. In both cases, each data point Pi is 

assigned a “theoretical point” or “model point” Pi

ˆ

 located on the proposed virtual object. That 

theoretical point is seen as the desired “correct” point, and the Euclidean distance between the 

two points is considered the individual error 

 i i ie = Pˆ − P 

of the data point with respect to the current location and shape specification of the virtual object. 

 

In orthogonal fitting, the theoretical point Pi

ˆ

 is chosen as a point that lies on the virtual object 

and is closest to the data point Pi in terms of Euclidean distance. In the 3D imaging 

environment, however, the data point Pi is considered to lie on a particular “scan ray” or “lineof-sight”, which emanates from the instrument position. 

In directional fitting, if the scan ray intersects the virtual object, the intersection closest to the

instrument is thus chosen as the theoretical point Pi

ˆ

 for the data point Pi. What happens if the 

scan ray of a data point Pi does not intersect the current virtual object? It might be tempting to 

reject such an occurrence as unrealistic as the point cloud was generated from the real object. It 

should be kept in mind, however, that during the fitting process, the virtual object will, in 

general, not match the actual object. Indeed, establishing that match is the purpose of the fitting 

process. It is, therefore, necessary to extend the error definition to those date points whose scan 

rays miss the virtual object. The following generic principle for a continuous extension has been 

proposed in [21]. Here, the theoretical point Pi

ˆ

 is chosen as a point on the virtual object that is

closest to the scan ray in terms of Euclidean distance. 

2.2 Sensitivity 

As we return to the general error function E (2.1.1), we examine a major aspect of analyzing the 

results of a fitting procedure. It concerns the “sensitivities” of the resulting parameters, namely, 

their marginal rates of change caused by perturbations of the data coordinates. Such sensitivities 

not only provide key information about a fitting process, they also play a role in the estimation of 

variances and covariances of the fitted parameters, as will be discussed in the subsequent section. 

With each set of stipulated data values d i n i i i,φ ,θ , ..., , , the error function E associates a set of 

minimizing parameters. We may thus consider the minimizing parameters as functions of these 

data variables 

 (2.2.1) ) ( , , ..., , , ,

1 1 1 n n n X d φ θ d φ θ
5

 ) ( , , ..., , , ,

1 1 1 n n n Y d φ θ d φ θ

 ) ( , , ..., , , ,

1 1 1 n n n Z d φ θ d φ θ

in a suitable neighborhood of the actually measured values d i n i i i, , , ...,1 ,

)0( )0( )0( φ θ = . By 

definition, 

 ) ( , , ..., , , ,

)0( )0( )0( )0(

1

)0(

1

)0(

1

)0(

n n n X = X d φ θ d φ θ

 ) ( , , ..., , , ,

)0( )0( )0( )0(

1

)0(

1

)0(

1

)0(

n n n Y = Y d φ θ d φ θ

 ) ( , , ..., , , ,

)0( )0( )0( )0(

1

)0(

1

)0(

1

)0(

n n n Z = Z d φ θ d φ θ

 are the desired results of the fitting process. 

In what follows, we will assume that the error function E satisfies all necessary differentiability 

conditions. We are particularly interested in the derivatives 

 , , , , , , , , , i ...,1 ,n,

X Y Z X Y Z

d

Z

d

Y

d

X

i i i i i i i i i

=

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

φ φ φ θ θ θ

because their values, for , , , , ...,1 ,

)0( )0( )0(

d d i n i = i φi = φi θ i = θ i = 

 (2.2.2) 

 

 ( ..., , ), ( ..., , ), ( ..., , )

)0( )0(

1

)0( )0( )0(

1

)0( )0( )0(

1

)0(

n

i i

n

i i

n

i i

d

X X

d

X X

d

d

X

d

X

θ

θ θ

θ

φ φ

θ

∂

∂

=

∂

∂

∂

∂

=

∂

∂

∂

∂

=

∂

∂

 ( ..., , ), ( ..., , ), ( ..., , )

)0( )0(

1

)0( )0( )0(

1

)0( )0( )0(

1

)0(

n

i i

n

i i

n

i i

d

Y Y

d

Y Y

d

d

Y

d

Y

θ

θ θ

θ

φ φ

θ

∂

∂

=

∂

∂

∂

∂

=

∂

∂

∂

∂

=

∂

∂

 ( ..., , ), ( ..., , ), ( ..., , )

)0( )0(

1

)0( )0( )0(

1

)0( )0( )0(

1

)0(

n

i i

n

i i

n

i i

d

Z Z

d

Z Z

d

d

Z

d

Z

θ

θ θ

θ

φ φ

θ

∂

∂

=

∂

∂

∂

∂

=

∂

∂

∂

∂

=

∂

∂

 

 

represent the respective sensitivities of the parameters )0( )0( )0( X , Y , Z to perturbations of the 

indicated data variables. 

 

Implicit differentiation will be used to derive expressions for the sensitivities (2.2.2) from the 

expression for the error function E . Indeed, the gradient of E with respect to these parameters, 

(2.2.3) 

































∂

∂

∂

∂

∂

∂

∇ =

( , , , , , ..., , , , )

( , , , , , ..., , , , )

( , , , , , ..., , , , )

1 1 1

1 1 1

1 1 1

n n n

n n n

n n n

XYZ

E X Y Z d d

Z

E X Y Z d d

Y

E X Y Z d d

X

E

φ θ φ θ

φ θ φ θ

φ θ φ θ
6

vanishes if the parameters X , have been minimized with respect to the coordinat Y, Z es 

d i n i i i,φ ,θ , = ...,1 , . As we thus substitute for the parameters X , their corresponding Y, Z

functions (2.2.1) in the above gradient components, we arrive at a set of derivative functions 

which are identically zero as functions of the data variables d i n i i i,φ ,θ , = ...,1 , . In other words, 

the following derivative expressions, 

 ( ( ..., , ), ( ..., , ), ( ..., , ), ..., , ) 0 1 1 1 1 ≡

∂

∂

n n n n X d Y d Z d d

X

E

θ θ θ θ 

(2.2.4) ( ( ..., , ), ( ..., , ), ( ..., , ), ..., , ) 0 1 1 1 1 ≡

∂

∂

n n n n X d Y d Z d d

Y

E

θ θ θ θ 

 ( ( ..., , ), ( ..., , ), ( ..., , ), ..., , ) 0 1 1 1 1 ≡

∂

∂

n n n n X d Y d Z d d

Z

E

θ θ θ θ 

vanish identically. Then so do the derivatives of these functions, which by the Chain Rule 

become: 

 

 0

2 2 2

2

2

≡

∂ ∂

∂

+

∂

∂

∂ ∂

∂

+

∂

∂

∂ ∂

∂

+

∂

∂

∂

∂

d X

E

d

Z

X Z

E

d

Y

X Y

E

d

X

X

E

i i i i

 0

2 2

2

2 2

≡

∂ ∂

∂

+

∂

∂

∂ ∂

∂

+

∂

∂

∂

∂

+

∂

∂

∂ ∂

∂

d Y

E

d

Z

Y Z

E

d

Y

Y

E

d

X

Y X

E

i i i i

 0

2

2

2 2 2

≡

∂ ∂

∂

+

∂

∂

∂

∂

+

∂

∂

∂ ∂

∂

+

∂

∂

∂ ∂

∂

d Z

E

d

Z

Z

E

d

Y

Z Y

E

d

X

Z X

E

i i i i

. 

For brevity, we displayed these relationships only for sensitivities with respect to the range 

variables id . Evaluating these functions for the resulting parameters 0 0 0 X ,Y ,Z and the actual 

data 0 0 0

, ,

i i i

d φ θ , yields a numerical m× m linear system of equations for the sensitivities (2.2.2) 

with respect to the coordinates id . With the notations 

( , , , ..., , ), ( , , , ..., , )

)0( )0(

1

)0( )0( )0(

2

)0(

2

)0( )0(

1

)0( )0( )0(

2

2

)0(

2

2

n n X Y Z d

X Y

E

X Y

E

X Y Z d

X

E

X

E

θ θ

∂ ∂

∂

=

∂ ∂

∂

∂

∂

=

∂

∂

, etc. 

this linear system takes the form 

 

)0(

2

)0( )0(

2

)0( )0(

2

)0( )0(

2

2

d X

E

d

Z

X Z

E

d

Y

X Y

E

d

X

X

E

i i i i

∂ ∂

∂

= −

∂

∂

∂ ∂

∂

+

∂

∂

∂ ∂

∂

+

∂

∂

∂

∂

(2.2.5) )0(

2

)0( )0(

2

)0( )0(

2

2

)0( )0(

2

d Y

E

d

Z

Y Z

E

d

Y

Y

E

d

X

Y X

E

i i i i

∂ ∂

∂

= −

∂

∂

∂ ∂

∂

+

∂

∂

∂

∂

+

∂

∂

∂ ∂

∂

 .

)0(

2

)0( )0(

2

2

)0( )0(

2

)0( )0(

2

d Z

E

d

Z

Z

E

d

Y

Z Y

E

d

X

Z X

E

i i i i

∂ ∂

∂

= −

∂

∂

∂

∂

+

∂

∂

∂ ∂

∂

+

∂

∂

∂ ∂

∂

 
7

The matrix of this linear system may be stated in terms of the Hessian 

(2.2.6) HXYZ E =

































∂

∂

∂ ∂

∂

∂ ∂

∂

∂ ∂

∂

∂

∂

∂ ∂

∂

∂ ∂

∂

∂ ∂

∂

∂

∂

2

2 2 2

2

2

2 2

2 2

2

2

Z

E

Z Y

E

Z X

E

Y Z

E

Y

E

Y X

E

X Z

E

X Y

E

X

E

of the error function E. The linear system (2.2.5) may thus be written as 

(2.2.7) )0(

)0(

)0(

)0(

)0( E

d

d

Z

d

Y

d

X

E XYZ

i

i

i

i

XYZ ∇

∂

∂

= −





































∂

∂

∂

∂

∂

∂

H , 

 

where again the symbol )0(is meant to indicate the, -- a posteriory --, substitution by 

)0( )0( )0( X , Y , Z and the actual data values. The resulting Hessian matrix is positive definite, and 

therefore nonsingular, at any locally unique minimum of the error function E. The linear system 

is then solvable and yields the values of the sensitivities (2.2.2) with respect to the variables i d . 

The remaining sensitivities with respect to the variables φi and θi may be determined from 

analogous linear systems, based on the same Hessian matrix. 

Note that implicit differentiation can be used to determine higher order sensitivities such as 

 

, , , ,, , ...,

2

2 2

2

2 2 2

2

2

i l j i j i l j l

X X X

d

X

d

X

d

X

φ θ φ φ θ ∂θ

∂

∂ ∂

∂

∂

∂

∂ ∂

∂

∂ ∂

∂

∂

∂

 . 

The corresponding linear systems are based on the same Hessian matrix as in (2.2.7) but use 

different right hand sides. 

2.3 Noise Propagation 

In, general, data variables i i id , and the parameters φ ,θ X , will be considered random Y, Z

variables with expected values )0(

1

)0( )0(

di,φi,θ and , , ,

)0( )0( )0( X Y Z respectively. In frequent 

applications, however, some data variables of the error function will be given “control variables” 

or “design variables”, and are thus not random. When fitting scanned objects, in particular, it is 
8

frequently assumed that the noise in range measurements id dominates noise in bearings, which 

furthermore is difficult to assess. Consequently, only the range coordinates id are considered 

random, while the bearing angles φi and θi are specified control variables. For scanning 

instruments, it is generally safe to assume that range variables id are independent of each other. 

The following exposition will be based on these assumptions, again for brevity. 

 

The sensitivities described in the previous section will be instrumental in assessing the effects of 

data noise on fitted parameters. The well known “Error Propagation Formula” provides first 

order estimates of the variances (see GUM [24], Chapter 5) 

 ∑=













∂

∂

≅

n

i

i

i

d

d

X

X

1

2

)0(

var( ) var( )

 (2.3.1) ∑=













∂

∂

≅

n

i

i

i

d

d

Y

Y

1

2

)0(

var( ) var( )

 ∑=













∂

∂

≅

n

i

i

i

d

d

Z

Z

1

2

)0(

var( . ) var( )

Similarly, one has for the covariances [25], 

 cov( , ) var( ) cov( , )

1

)0( )0(

d Y X

d

Y

d

X

X Y i

n

i i i

 =











∂

∂

∂

∂

≅ ∑=

 (2.3.2) cov( , ) var( ) cov( , )

1

)0( )0(

d Z Y

d

Z

d

Y

Y Z i

n

i i i

 =











∂

∂

∂

∂

≅ ∑=

 cov( , ) var( ) cov( , )

1

)0( )0(

d X Z

d

X

d

Z

Z X i

n

i i i

 =











∂

∂

∂

∂

≅ ∑= . 

Again, it is here assumed that the coordinate measurements )0(

i

d are independent, that is, not 

correlated, so that their mutual covariances are zero. Note however that, even if the measured 

quantities )0(

i

d are independent, the fitted parameters )0( )0( )0( X , Y , Z will still be correlated. 

In most applications, the condition of homoscedacity is supposed to hold, that is, the variances 

have the same value within a class of measurements, such as the class d of range measurements, 

 d d i n i

var( ) = var( ), = ...,1 . 

For the special cases of LS and NLS regression, in which the individual errors in our scenario 

would take the form 

 ) ( , , , ,

i i i X Y Z i i

e = d − f φ θ , 
9

a covariance matrix can be approximated in an elegantly simple fashion under homoscedacity as 

set forth in the general literature, e.g. [5-7, 26, 27]. Unfortunately, the error functions considered 

here, in particular, the orthogonal error function (Chapter 4) and a portion of the directional error 

function (Chapter 3) do not fall into this regression category as the required separation of the 

random variables from the control portion cannot be achieved. It is for this reason, that the more 

general approach described in Sections 2.2 and 2.3 had to be adopted. A more detailed 

explanation will be provided in a separate report [23]. 

3 Directional Fitting of Spheres 

3.1 Directional Errors 

Introducing the trigonometric quantities i i i ξ , , the Cartesian coordinates η ,ς i i ix , of data y ,z

points will be expressed in the form 

(3.1.1) i i i i i i

x = d cos , φ cosθ = d ξ i i i i i iy = d sin , φ cosθ = d η i i i i iz = d sin , θ = d ς

where 1

2 2 2

ξ i + ηi + ς i = . The vector ) ( , ,

i i i ξ η ς represents the direction of the scan ray along 

which the data point Pi was acquired. Next we introduce the quantities: 

 

 i Xi i Yi i Zi ip = ξ + η + ς , 

(3.1.2) 2 2 2 2 2qi = X + Y + Z − pi, 

 

2 2 2

i R qi

s = − . 

Figures 1 and 2 illustrate the geometric meaning of these quantities. 
10

If the scan ray of data point Pi

 intersects the virtual sphere centered at ] C = [X Y Z , we 

associate with Pi

 the midpoint Pi

′ of the resulting chord. The quantity pi Pi = ′ then 

represents the distance of that “mid-chord” point from the instrument at the origin ]0 O = 0[ 0 . 

Similarly, qi = P′ − C ≥ 0 represents the distance of the mid-chord from the sphere center. 

Thus 

 qi < R

is the condition for true, that is, non-tangential intersection. The quantities i

p and iq are side 

lengths of the right triangle ∆OPi

′C , with its right angle at Pi′ . Pythagoras thus yields the 

relation (3.1.2) between 2

i

p and 2

i

q . The triangle ∆ Pˆ

i Pi

′C , where the theoretical point Pi

ˆ

 

marks the entry point into the sphere, also has a right angle at Pi

′ . The quantity i

s thus 

represents the length of the half-chord that needs to be subtracted from the distance = P′ pi to 

arrive at the desired distance Pˆ of the theoretical point from the origin. As a result, the 

directional error of the data point Pi

 is given by 

(3.1.3) i i i i

f = p − s − d (= “interior error” of Pi if qi < R ) . 

pi

di

qi

C

R

si

Pi

Pi

^

O

Pi

p ’

i

di

qi

C

R

si

Pi

Pi

^

O

Pi

’

Figure 1. Geometrical interpretation of the directional error function when a scan ray 

intersects the sphere surface. Pi

 (marked by dark dot) is the experimental point, light dots 

mark the theoretical point on the sphere surface Pi

ˆ

 and the mid-chord point Pi

′ . The length 

of the bold line segment measures the error if defined by (3.1.3). 
11

On the other hand, if 

 qi ≥ R , 

then the scan ray of the data point Pi fails to truly intersect the virtual sphere. In that case, and 

following the general extension principle set forth in Section 2.1, we determine the theoretical 

point as that point on the virtual sphere which is closest to the scan ray. The line segment which 

represents the shortest distance between the sphere and the scan ray has to be orthogonal to both 

the sphere and the scan ray. The line segment thus has to be part of a line through the center of 

the sphere, and also meet the scan ray at a point Pi

′ at a right angle. This defines again the right 

triangle ∆OPi

′C , which we encountered before, and whose side lengths are again i

p and iq

(Fig.2). The desired theoretical point Pi

ˆ

 thus lies on the side [ Pi

′O ] at distance R from center 

and at distance qi − R from Pi

′ . The triangle ∆ Pˆ

i Pi

′C is also a right triangle, has side lengths 

i i i i P ′ − P = p − d and Pi

′− Pˆ

i = qi − R . The length of its hypotenuse Pi − Pi

ˆ

 thus 

represents the error of the data point Pi : 

(3.1.4) 2 2

g ( p d ) (q R) i = i − i + i − (= ”exterior error” of Pi if qi ≥ R ). 

pi

di qi

C

Pi

Pi

^

O

Pi

’

pi

di qi

C

Pi

Pi

^

O

Pi

’

Figure 2. Geometrical interpretation of the directional error function when a scan ray does 

not intersect the sphere surface. Pi (marked by dark dot) is the experimental point, light dots 

mark the theoretical point on the sphere surface Pi

ˆ

 and the point Pi

′ on a scan ray which is 

the closest one to the sphere center C. The length of the bold line segment measures the

error ig defined by (3.1.4). 
12

(Note that, by Thales’ theorem, the locus of all possible points Pi

′ is the sphere through both the 

virtual sphere center C and the originO , centered halfway between these two points.) 

If 0 qi − R = then i if = g , so that the combined error function will be continuous. While the 

error expression ig is everywhere twice continuously differentiable, the error expression if fails 

to be so if and only if 0 si = , -- the case of tangential intersection -- , where its gradient with 

respect to the parameters X , is infinite. In these cases, the resulting full Y, Z error function 

will also not be differentiable. However, those points will only amount to a closed set of 

measure zero in parameter space. As a consequence, a gradient based numerical minimization 

method, such as the often relied on “BFGS” method [28], may still be used [21]. Similarly, the 

probability of the error function E not being differentiable for the fitted parameters

)0( )0( )0( X , Y , Z will be theoretically zero. 

We find it convenient, to categorize only a true intersections as a “hit”. A tangential intersection 

is thus considered a “miss”, along with all cases in which the virtual sphere is not met at all. 

Accordingly, we divide the indices i into two sets: 

 U = {i : and qi < R} V = { j : . qj ≥ R}

The combined error function then takes the form 

 (3.1.5) = = ∑i V ∈ ∈ U+ ∑j

direc direc n n n i g j E E X Y Z d d f

2 2

1 1

( , , , ,φ ..., , , φ , θ )

3.2 Derivatives for the Directional Error Function

In this section, we list formulas for the gradients and the Hessians of the individual error 

functions i

f and ig with respect to the parameters X , , along with the second derivatives Y, Z

with respect to both these parameters and the data variables i i id , . Gradients support φ , θ

optimization methods and are the first step towards determining the above second derivatives, 

which are needed for the computation of the sensitivities and variances described in Sections 2.3-

4. Derivation of these formulas is provided in the Appendix as referenced. 

In terms of the individual errors if and ig , the gradients and Hessians of the directional error 

function are: 

 ∇ = ∑∈ ∈∇ + ∑ ∇

i V U i XYZ drec XYZ i XYZ i E f g

2 2

 = ∑i V ∈ ∈ U

+ ∑i XYZ drec XYZ i XYZ i E f g

2 2 H H H . 

The gradients are linear combinations of the vectors 
13

(3.2.1) 





























=

Z

Y

X

U and 





























=

i

i

i

i

ς

η

ξ

Λ , 

namely (see (A.2.6) and (A.2.7) in Appendix) 

(3.2.2) i

i

i i i

i

i

XYZ i

s

f f d

s

f

f U Λ

( )

2

1 2 +

∇ = − , i i

i

i

i

XYZ i d

q

Rp

q

R

g 1( )U ( )Λ

2

1 2 ∇ = − + −

 

Similarly, the Hessian matrices are linear combinations of the following four symmetric 

matrices: 

(3.2.3) 





























=

0 0 1

0 1 0

1 0 0

I (Identity), 

 [ ]





























=





























=

2

2

2

ZX ZY Z

YX Y YZ

X XY XZ

X Y Z

Z

Y

X

T UU , 

 

[ ] [ ]





























+ + +

+ + +

+ + +

=





























+





























+ =

Z X Z Y Z Z

Y X Y Y Y Z

X X X Y X Z

X Y Z

Z

Y

X

i i i i i i

i i i i i i

i i i i i i

i

i

i

i i

T

i

T

i

ξ ς η ς ς ς

ξ η η η ς η

ξ ξ η ξ ς ξ

ς

η

ξ

UΛ ΛU ξ η ς , 

 [ ]





























=





























=

2

2

2

i i i i i

i i i i i

i i i i i

i i

i

i

i

T

i i

ς ξ ς η ς

η ξ η η ς

ξ η ξ ξ ς

ξ η ς

ς

η

ξ

Λ Λ . 

Thus, (see (A.2.11) and (A.2.12)) 
14

(3.2.4) 

T

i i

i

i i

i

i i i i

T

i

T

i

i i

T i i i

i

i i

i

i i

XYZ i

s

p s

s

p s p d

s s

p p d

s

p d

s

p d

f

Λ Λ

H I UU UΛ ΛU















 −

+

− −

+

+ 















−

−

− 













 −

+ 















−

−

=

( )( ) (2 )

( )

( ) 1

1

2

1

3

2 2

3 3

2

 

T

i i

i

T

i

T

i

i

T i

i i

XYZ i

q

X Y Z

R

q

Rp

q

R

q

R

g

Λ Λ

H I UU UΛ ΛU















 + +

+

+ − + 















= −

3

2 2 2

3 3

2

1 ( )

2

1

. 

Evaluated for optimal parameters )0( )0( )0( X , Y , Z and actual data points ] [

)0( )0( )0(

i i i i P = d φ θ , 

these Hessian matrices support the left hand side of linear systems (2.2.7) for the corresponding 

sensitivities. For the right-hand sides of those systems, we have for the range variables i d , 

 ∑ ∑ ∈ ∈ ∂

∂

+ ∇

∂

∂

= ∇

∂

∂

∇

i V U i

i

i

XYZ

i

i

drec XYZ

i

XYZ d

g

d

f

E

d

2 2

where (see (A.3.1) and (A.3.2)) 

(3.2.5) i

i

i i

i i

i

XYZ

s

p s

d s

f

U Λ

−

= − −

∂

∂

∇

1

2

1

2

 , i

i

i

XYZ d

g

= − Λ

∂

∂

∇

2

2

1

 . 

For the bearing variables φ , θ , the individual derivatives in 

 ∑ ∑ ∈ ∈ ∂

∂

+ ∇

∂

∂

= ∇

∂

∂

∇

i V U i

i

i

XYZ

i

i

drec XYZ

i

XYZ

f g

E

φ φ φ

2 2

 , 

 ∑ ∑ ∈ ∈ ∂

∂

+ ∇

∂

∂

= ∇

∂

∂

∇

i V U i

i

i

XYZ

i

i

drec XYZ

i

XYZ

f g

E

θ θ θ

2 2

. 

are multiples of the vectors 

(3.2.6) 



























−

=

∂

∂

0

i

i

i

i

ξ

η

φ

Λ , 





























=

∂

∂

i

i

i

i

i

γ

β

α

θ

Λ , 

where i i i i i i i θi α = −cosφ sinθ , β = −sinφ sinθ , γ = cos . Thus 
15

(3.2.7) 



























−

=

∂

∂

∇

0

( )

2

1 2

2

i

i

i i

i

i

XYZ f

f

ξ

η

φ

Γ , 



























−

=

∂

∂

∇

0

( )

2

1 2

2

i

i

i i

i

i

XYZ g

g

ξ

η

φ

Γ

(3.2.8) 





























=

∂

∂

∇

i

i

i

i i

i

i

XYZ f

f

γ

β

α

θ

( )

2

1 2

2

Γ , 





























=

∂

∂

∇

i

i

i

i i

i

i

XYZ g

g

γ

β

α

θ

( )

2

1 2

2

Γ , 

where the two pre-multiplying matrices are given by

(3.2.9) 

( ) 1 ,

( )( ) 2

( )

1 ( )

( )

3 2

2 Γ UU ΛU I















 −

+ − − 















−

+ −

+ − 













 −

= −

i

i i

i i

T

i

i i

i i i i

i i

T

i

i i i

i

i i

s

p d

p s

s s

p s p d

p s

s

p p d

s

f

 Γ UU ΛU I

















+ − 













 + +

+ 















= − u

i

T i

i

i

T

i

i

i i d

q

Rp

q

X Y Z

R

q

Rp

g 3

2 2 2

3

2

( ) . 

As 

(3.2.10) Λ = 0



























−

0

i

i

T

i ξ

η

 and Λ = 0





























i

i

i

T

i

γ

β

α

 , 

the above matrices may be replaced in (3.2.7) and (3.2.8), respectively, by the following 

symmetrized versions: 

(3.2.11) 

( ) 1 ,

( )

( )( ) 2

( )

1 ( )

( )

3 2

2

I

Γ UU ΛU UΛ















 −

− −

+ + 















−

+ −

+ − 













 −

= −

i

i i

i i

T

i

T

i

i i

i i i i

i i

T

i

i i i

i

i

S

i

s

p d

p s

s s

p s p d

p s

s

p p d

s

f

 
16

 Γ UU ΛU UΛ I

















+ + − 













 + +

+ 















= − u

i

T i

i

T

i

i

T

i

i

i

S

i d

q

Rp

q

X Y Z

R

q

Rp (g ) ( )

3

2 2 2

3

2

. 

4 Orthogonal Fitting of Spheres 

Here, the theoretical point Pi

ˆ

 , that is, the point on the sphere which is closest to the data point 

[ ] i i i i P = x y z , defines the individual error 

 hi = Pˆ

i − Pi = C − Pi − R = wi − R

with respect to the sphere center C and the radius R . We thus represent the orthogonal fitting 

approach by the full error function 

(4.1.1) 2 2

1

2

1

2 E h (w R) wi 2Rwi R

n

i

i

n

i

orth = ∑ I = ∑ − = − +

= =

 . 

Consistent with the generation of point clouds by scanning from a single instrument location, and 

as discussed before, the underlying coordinate frames are again considered polar with the 

instrument location at the origin. For an analytic discussion of the orthogonal error function in 

terms of Cartesian data see [11]. 

With the notation (3.1.1) and the definition (3.1.2) of the auxiliary quantitity i p , 

 i i i i i i

x = d cos , φ cosθ = d ξ i i i i i iy = d sin , φ cosθ = d η i i i i iz = d sin , θ = d ξ

we have 

(4.1.2) 2 2 2 2 2 2 2 2( ) ( ) ( ) ( ) 2 i i i i i i i w = X − x + Y − y + Z − z = X + Y + Z − d p + d . 

A key vector, in which gradients and Hessians of the individual orthogonal error squares 2

i h may 

be expressed, is given by 

(4.1.3) i i

i

i

i

i

i

i

i

i d d

Z

Y

X

Z z

Y y

X x

W = U − Λ





























−





























=





























−

−

−

=

ς

η

ξ

, 
17

as 2

i i

T Wi W = w . Also (see (A.4.3) and (A.4.5)), 

(4.1.4) i

i

XYZ i

w

R

h



W













∇ = 1−

2

1 2

and 

(4.1.5) T

i i

i i

XYZ i WW

w

R

w

R

h 3

2

1

2

1

+ 















H = − I . 

For the derivatives which define the right-hand sides of the linear system (2.2.7), we have first: 

 ∑= ∂

∂

= ∇

∂

∂

∇

n

i i

i

orth XYZ

i

XYZ d

h

E

d 1

2

with (see (A.5.2)) 

(4.1.6) i

i

i i i

i i

i

XYZ

w

R

p d

w

R

d

h W Λ















= − − − −

∂

∂

∇ ( ) 1

2

1

3

2

 . 

Again, the corresponding mixed derivatives with respect to the bearing variables φi θi

, 

are multiples of the vectors i

i

Λ

∂φ

∂

and i

i

Λ

∂θ

∂

, defined in (3.2.6). Their common multiplier is the 

matrix 

 (4.1.7) T

i

i

i

i

i i i

w

Rd

w

R

Γ h d I

3 W U

2

( ) 1 + 















= − − , 

which, in analogy to (3.2.11), may be replaced also by its symmetrized form, 

(4.1.8) T

i i

i

i

i

i i

S

i

w

Rd

w

R

Γ h d I

3 WW

2

( ) 1

2

1

+ 















= − − . 

Thus 

(4.2.5) 



























−

=



























−

=

∂

∂

∇

0

( )

0

( )

2

1 2 2

2

i

i

i

S

i i

i

i i

i

i

XYZ h h

h

ξ

η

ξ

η

φ

Γ Γ

and 
18

(4.2.6) 





























=





























=

∂

∂

∇

i

i

i

i

s

i

i

i

i

i i

i

i

XYZ h h

h

γ

β

α

γ

β

α

θ

( ) ( )

2

1 2 2

2

Γ Γ . 

This concludes the main part of the report. It is followed by the Appendix in which details about 

the derivation of the key formulas are provided. 

Appendix A: Determination of Derivative Formulas Used for Calculating Sensitivities 

Here, we provide step by step developments of the derivative formulas referred to in Chapters 3 

and 4 for the purpose of determining the parameter sensitivities for directional and orthogonal 

fitting. In Section A.2, the gradient ∇XYZEdirec and Hessian HXYZEdrec of the directional error 

function are at issue. The Hessian provides the matrix for the corresponding linear system 

(2.2.7). Also for the directional error function, the derivatives of both parameters and data 

variables, 

 

i

XYZ XYZ

i i

XYZ XYZ

i i

XYZ XYZ

i

E

E

E

E

d

E

E

d φ φ θ ∂θ

∂

∇ = ∇

∂

∂

∂

∂

∇ = ∇

∂

∂

∂

∂

∇ = ∇

∂

∂

, , , 

are derived in Section A.3, furnishing the right hand sides of these systems. Finally, Section A.4 

provides the analogous information in the case of orthogonal fitting. 

A.1 General Considerations 

In what follows, the calculation of gradients ∇XYZ and Hessians HXYZ will often be based on the 

following straightforward reformulations of product and chain rules: 

(A.1.1) ∇XYZb(a) = b′(a)∇XYZ a , 

 ∇XYZ ab = a∇XYZb + b∇XYZ a , 

 a a a ∇XYZ = ∇XYZ 2

2

 , 

 

2 1 2

2

1

a

a

a a ∇XYZ = ∇XYZ = ∇XYZ , 

and 

(A.1.2) b a b a a a b a a XYZ

T HXYZ ( ) = ′′( )∇XYZ ∇XYZ + ′( )H , 

 ab a b b a a b b a XYZ XYZ

T

XYZ XYZ

T HXYZ = ∇XYZ ∇XYZ + ∇ ∇ + H + H , 
19

 a a a a a XYZ

T HXYZ 2 XYZ XYZ 2 H2

= ∇ ∇ + , 

 

2 2 2

3

2 1

2

1 1

4

1

a

a

a a

a

a a XYZ

T HXYZ = HXYZ = − ∇XYZ ∇XYZ + H . 

These formulas are straightforward reformulations of product and chain rules, and will not 

always be referred in what follows. 

A.2 Gradients and Hessians of the Directional Error Functions 

Recall the directional error function 

 = = ∑i V ∈ ∈ U

+ ∑i

drec i i E f g

2 2

with individual errors (3.1.3-4), 

 i i i if = p − s − d , 

2 2

g ( p d ) (q R)

i = i − i + i − , 

based on the auxiliary quantities i i i p , (3.1.2) and the direction cosines (3.1.1) q , s

 i i i i i i i θi ξ = cosφ cosθ , η = sinφ cosθ , ς = sin . 

Again, all gradients ∇XYZ determined in this Section will be linear combinations of the two 

vectors, 

 





























=

Z

Y

X

U and 





























=

i

i

i

i

ς

η

ξ

Λ , 

which were introduced in Section 3.2. Using (A.1.1) where indicated, we note: 

(A.2.1) XYZ i XYZ X i Y i Z i i ∇ p = ∇ ( ξ + η + ς ) = Λ , XYZ pi 2 piΛi

2 ∇ = , 

 ( ) [2 ]

2 2 2 2 2

XYZ i XYZ XYZ i i i ∇ q = ∇ X + Y + Z − ∇ p = U − p Λ , 

 [ ] i i

i

XYZ i

i

XYZ i XYZ i

p

q

q

q

∇ q = ∇ q = ∇ = U − Λ

1 1

2

2 1 2

 , 
20

 [ ] XYZ i XYZ i XYZ i i i ∇ s = ∇ (R − q ) = −∇ q = −2 U − p Λ

2 2 2 2

 , 

 [ ] i i

i

XYZ i

i

XYZ i XYZ i

p

s

s

s

∇ s = ∇ s = ∇ = − U − Λ

1 1

2

2 1 2

 , 

 [ ] i i

i

XYZ i

i i

XYZ p

s

s

s s

∇ = − 2 ∇ = 3 U − Λ

1 1 1

 . 

All Hessians HXYZ determined in this section will be linear combinations of the four matrices 

(3.2.4). Again we begin with the auxiliary quantities: 

 

(A.2.2) 0 HXYZ pi = , 

T

i i i

T ΗXYZ pi 2 XYZ pi XYZ p 2 Λ Λ2

= ∇ ∇ = by (A.1.2), 

 ( )

T

XYZ i XYZ XYZ i i i H q = H (X + Y + Z ) − H p = 2 I − Λ Λ

2 2 2 2 2

 , 

 

 ( )

T

XYZ i XYZ i i i H s = − H q = −2 I − Λ Λ

2 2

 , 

From Hessians of 2

qi

, 

2

i

s we pass to Hessians of i

q , is , using (A.1.2) 

 ( )( ) ( )

( .)

1

( )

1

1 1

1

2

1 1

4

1

3

2

3 3

3

2 2 2

3

2

T

i i

i

T

i i

i

i

i i

i

T i

i

T

i i

i

T

i i i i

i

XYZ i

i

i

T

XYZ i XYZ

i

XYZ i XYZ i

q q

p

q

p

q

q

p p

q

q

q

q q

q

q q

UU UΛ ΛU Λ Λ I Λ Λ

U Λ U Λ I Λ Λ

H H H

= − + + − + −

= − − − + −

∇ ∇ +

−

= =

Thus 

(A.2.3) T

i i

i i

i

i i

i

T i

i i

XYZ i

q q

p

q

p

q q

H q I UU UΛ ΛU Λ Λ 















= − + + − +

1

( )

1 1

3

2

3 3

. 

Concerning the last term, note 

 3

2 2 2

3

2 2

3

2

1

i i

i i

i i

i

q

X Y Z

q

p q

q q

p + +

=

+

+ = . 

Similarly, by (A.1.2), 
21

 ( )( )

( ). 1( )

1

( )

1 1

1

2

1 1

4

1

3

2

3 3

3

2 2 2

3

2

T

i i

i

T

i i

i

i

i i

i

T i

i

T

i i

i

T

i i i i

i

XYZ i

i

i

T

XYZ i XYZ

i

XYZ i XYZ i

s s

p

s

p

s

s

p p

s

s

s

s s

s

s s

UU UΛ ΛU Λ Λ I Λ Λ

U Λ U Λ I Λ Λ

H H H

= − + + − − −

= − − − − −

∇ ∇ +

−

= =

Thus 

(A.2.4) T

i i

i

i

i

i i

i

T i

i i

XYZ i

s

p

s s

p

s s

H s I UU UΛ ΛU Λ Λ 















= − − + + + − 3

2

3 3

1

( )

1 1

. 

Concerning the last term, note 

 3

2 2 2 2

3

2 2 2 2

3

2 2 2

3

2 2

3

2

1 ( )

i i i

i i

i

i i

i

i

i

q

X Y Z R

q

R X Y Z

q

R q p

q

s p

s

p

s

+ + −

= −

− + +

=

− −

=

−

− = . 

From the above, we derive derivative expressions involving the errors if and ig : 

 

(A.2.5) 

( )

i

i

i i

i

i

i

i i

i

i i

i

XYZ i XYZ i XYZ i i

s

f d

s s

p s

s

p

s

f p s

U Λ U Λ

Λ U Λ

+

= −

−

= −

∇ = ∇ − ∇ = + −

1 1

1

(A.2.6) i

i

i i i

i

i

XYZ i i XYZ i

s

f f d

s

f

f f f U Λ

( )

2

1 2 +

∇ = ∇ = −

For the external portion of the error function, we find 

( )

i

i

i i i i i

i

i

i i

i

i

XYZ i i i XYZ i i XYZ i i i

q

p d q q R p

q

q R

p

q

q R

g p d p q R q p d

U Λ

U U Λ

(2 ) (2 ) (2 )

(2 )

(2 ) (2 ) (2 )

2

− − −

+

−

=

−

−

∇ = − ∇ + − ∇ = − +

or 

(A.2.7) i i

i

i

i

XYZ i d

q

Rp

q

R

g 1( )U ( )Λ

2

1 2 ∇ = − + −
22

and by (A.1.1), 

(A.2.8) ( ) i i i i i

i i

XYZ i

i

XYZ i XYZ i q R R p d q

q g

g

g

g g ( )U ( )Λ

1 1

2

2 1 2 ∇ = ∇ = ∇ = − + − . 

Moving to the second derivatives, we find 

(A.2.9) T

i i

i

i

i

T

i

T

i

i

T i

i i

XYZ i XYZ i

s

p

s s

p

UU

s s

H f H s I UΛ ΛU Λ Λ 















= − = + − + + − 3

2

3 3

1

( )

1 1

 . 

Next, we introduce the matrix 

 

( )( )

T

i i

i

T i i

i

T

i

i

T i i

i

T

i i i i i i

i

i

T

XYZ i XYZ i XYZ

s

p s

s

p s

s

p s p s

s

f f f

UU UΛ ΛU Λ Λ

G U Λ U Λ

2

2

2 2

2

( )

( )

1

( ) ( )

1

−

+ +

−

= −

= ∇ ∇ = − − − −

. 

Also by (A2.9), 

 

T

i i

i

T i i i

i

T

i

i

T i i

i

i

i

i

i XYZ i

s

f p s

s

f p

s

f

s

f

f H f I UU UΛ ΛU 3 Λ Λ

2 2

3 3

( )

( )

−

= + − + + . 

By (A.1.2) , XYZ i XYZ i i XYZ i H f 2G f 2 f H f

2

= + . Thus 

(A.2.10) 

.

( ) ( )

( )

1

2

1

3

2 2

2

2

2 3 2 3

2

T

i i

i

i i i

i

i i

T

i

T

i

i

i i

i

T i i

i

i

i i

i

XYZ i

s

f p s

s

p s

s

f p

s

p s

s

f

s s

f

f

Λ Λ

H I UU UΛ ΛU















 −

+

−

+

+ 















+

−

− 















= + +

 

Expressing i i i i

f = p − s − d yields 

2 3 3 3

1

i

i i

i

i i

i

i

i

s

p d

s

s f

s

f

s

−

=

+

+ =

i i

i i i

i

i i i i i i i i

i

i i i i i

i

i i

i

i i

s s

p p d

s

p s s p p s p d

s

p s s f p

s

f p

s

p s ( ) ( ) 1

3 3

2 2

2 3 3

−

−

=

− + − −

=

− +

+ =

−
23

3

2 2 2

3

3 2 2 2 2 3

3

2 2 3 3 2 2 3 2 2

3

2 2

3

2

( )( ) (2 )

( ) ( ) 2 2

( ) ( ) 2

i

i i i i i i i

i

i i i i i i i i i i

i

i i i i i i i i i i i i i i i

i

i i i

i

i i

s

p s p d p s s

s

p p d p s d s p s s

s

p s p s s p p s s p s p d d s

s

f p s

s

p s

− − − −

=

− − − − +

=

− + + − − + − +

=

−

+

−

From these relations, we find the alternate expression 

(A.2.11) 

.

( )( ) (2 )

( )

( ) 1

1

2

1

3

2 2

3

3

2

T

i i

i

i i

i

i i i i

T

i

T

i

i i

i i i

T

i

i i

i

i i

XYZ i

s

p s

s

p s p d

s s

p p d

s

p d

s

p d

f

Λ Λ

UΛ ΛU

H I UU















 −

+

− −

+

+ 















−

−

−















 −

+ 















−

−

=

which duplicates the first portion of (3.2.1). 

We move now to the external portion i

g of the error function. Note that the following Hessians 

vanish: 

 

2 2 HXYZ pidi

, HXYZ d i, HXYZ R . 

Thus 

( )

XYZ i i i XYZ i XYZ i

XYZ i XYZ i i i XYZ i i i i i i

p q Rq X Y Z Rq R q

g p d q R p p d d q q R R

H H I H

H H H

( 2 ) ( 2 ) 2 2

( ) ( ) ( 2 2

2 2 2 2 2

2 2 2 2 2 2 2

= + − = + + − = −

= − + − = − + + − +

and by (A.2.2), 

(A.2.12) T

i i

i

T

i

T

i

i

T i

i i

XYZ i

q

X Y Z

R

q

Rp

q

R

q

R

H g I UU UΛ ΛU Λ Λ 













 + +

+ − + + 















= − 3

2 2 2

3 3

2

1 ( ) . 

This establishes the second portion of (3.2.1). 

A.3 Mixed Derivatives of the Directional Error Function 

The right-and-sides of the system of linear equations (2.2.7) are at issue. They require the 

negatives mixed derivatives of the form 
24

 XYZEdirec XYZ Edirec ∂ ∗

∂

∇ = ∇

∂ ∗

∂

 , 

where ∗ indicates a data variable of the error function. We first consider the data variable i

d . 

Note, in this context, 

= ∇ (− )1 = 0

∂

∂

∇ XYZ

i

i

XYZ d

f

, ( ) XYZ i XYZ i XYZ i

i

i

XYZ i

i

i

XYZ f p s

d

f

f

d

f

= ∇ − = − ∇ − ∇ 















∂

∂

= ∇

∂

∂

∇ 2 ( 2 ) 2

2

 . 

Thus by (A.2.1), 

(A.3.1) ( ( ) )

1

( )

1

2

1

2

i i i

i

i i

i

i

i

i

XYZ p s

s

p

d s

f

= −Λ − U − Λ = − U + − Λ

∂

∂

∇ . 

Similarly, 

( ) ( ) i i i XYZ i i XYZ i

i

XYZ

i

i

XYZ p d q R p d p

d d

g

− + − = − ∇ − = − ∇

∂

∂

= ∇

∂

∂

∇ ( ) ( ) 2 ( 2

2 2

2

so that 

(A.3.2) i

i

i

XYZ d

g

= −Λ

∂

∂

∇

2

2

1

 . 

As pointed out in Section 3.2, the calculation of the corresponding derivatives with respect to φi

and θi will be based on the matrices generated by the differential operator 

 [ ] XYZ

T

i

i i i

Γ = ∇ξ η ς ∇

introduced in Section 3.2. In order to apply this operator to the individual errors 2 2,

i i

f g , we 

first apply the transposed gradient T

i i i

∇ξ η ς to auxiliary quantities. In particular, note: 

(A.3.3) 

 
25

Λ = I





























=





































∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∇ =

0 0 1

0 1 0

1 0 0

i

i

i

i

i

i

i

i

i

i

i

i

i

i

i

i

i

i

i

T

i i i

ς

ς

η

ς

ξ

ς

ς

η

η

η

ξ

η

ς

ξ

η

ξ

ξ

ξ

ξ η ς

, U = 0





























=





































∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∂

∇ =

0 0 0

0 0 0

0 0 0

i i i

i i i

i i i

T

Z Z Z

Y Y Y

X X X

i i i

ξ η ς

ξ η ς

ξ η ς

ξ η ς

and 

(A.3.4) [ ]

T

i i i

T

i

T

p X Y Z X Y Z i i i i i i∇ξη ς = ∇ξη ς( ξ + η + ς ) = = U , 

 

T

i i

T

i i

T

p p p p

i i i i i i

2 ( ) 2 U

2 ∇ξ η ς = ∇ξ η ς =

 

T

i i

T

i

T

i

T

q X Y Z p p p

i i i i i i i i i

( ) ( ) 2 U

2 2 2 2 2 2 ∇ξ η ς = ∇ξ η ς + + − = −∇ξ η ς = −

 

T

i

i

i

T

i

i

T

i

T

q

p

q

q

q q

i i i i i i i i i

∇ = ∇ = ∇ ( ) = − U

1

2

1

)

2

2

2

ξ η ς ξ η ς ξ η ς

 

T

i

i

i

T

i i

T

q

p

q

q q

i i i i i i 2

( )

3 U

1 1

∇ξ η ς = − ∇ξ η ς =

 

T

i i

T

i

T

i

T

s R q q p

i i i i i i i i i

( ) ( ) 2 U

2 2 2 2 ∇ξ η ς = ∇ξ η ς − = −∇ξ η ς =

 

T

i

i

i

T

i

i

T

i

T

s

p

s

s

s s

i i i i i i i i i

∇ = ∇ = ∇ ( ) = U

1

2

1 2

2

2

ξ η ς ξ η ς ξ η ς

 

T

i

i

i

T

i i

T

s

p

s

s s

i i i i i i 2

)

3 U

1 1

∇ξ η ς = − ∇ξ η ς = −

Consequently, 

 (A.3.5) T

i

T i i

i

i

i

T

i

T

i i i

T

i

T

s

f d

s

p

f p s d p s

i i i i i i i i i i i i

U U













 +

= − 















∇ξ η ς = ∇ξ η ς( − − ) = ∇ξ η ς − ∇ξ η ς = 1−

With these intermediary results, we are able to determine the matrices ) ( ,) (

2 2

i i i i Γ f Γ g . By 

(A.2.6) and the Product Rule, 
26

 

.

( ) ( )

( )

( )

2

1

( )

2

1 2 2















 +

− ∇

+

− ∇ 















= ∇















 +

= ∇ ∇ = ∇ −

i

T i i i

i

i

i i i

i

T

i

T i

i

i i i

i

i

T i

XYZ i

T

i i

s

f f d

s

f f d

s

f

s

f f d

s

f

f f

i i i i i i i i i

i i i i i i

ξ η ς ξ η ς ξ η ς

ξ η ς ξ η ς

U Λ Λ

Γ U Λ

 . 

By (A.3.3-5), 

T

i

T i i i i i i

i

i i i i i i i i

T

i

T i i i i i i i

i

i i

i

i i

T

i

i

i

i

T

i

i i

i

T

i

i

i

T

i

T i

s

f d f s sf

s

f d s f f d sf

s

f d s f f s d

s

f p

s

f d

s

p

f

s s

f d

s

f

s

f

s

f

i i i i i i i i i

U U

U U

U U















 + + +

= −















 + + + +

= −















 + + + +

= − 















−

+

= −

















+ −















 +

= ∇ + ∇ = −

















∇

3 3

2 3 3

3

( ) ( ) ( )( )

( ) ( )

1 1 1

( )

ξ η ς ξ η ς ξ η ς

( )

( ) T

i

i i i i i i i i

T

i i i i i i i i

i

i i

T

i

i i

i

i

i i

i

i i i i i i

i i

T

i

i

i i

i

T i

i

T i i i

s

f d f d f s sf

f d f s sf sf

s

f d

s

f d

s

f

f d

s

f d f s sf

f d

s

f

f d

s

f

s

f f d

i i i i i i i i i

U

U

U















 + + + +

= −

+ + + + 













 +

= −































 +

+ + 













 + + +

= −

+ + ∇ +

































= ∇ 













 +

∇

3

3

3

( () )( ) 2 )

( )( ) )

( )

( )( )

( ) ( )

( )

ξ η ς ξ η ς ξ η ς

From this and by (A.3.3), 

(A.3.6) 

( )

I

Γ UU ΛU















 +

−















 + + + +

− 













 + + +

= −

i

i i i

T

i

i

T i i i i i i i i

i

i i i i i i

i i

s

f f d

s

f d f d f s sf

s

f d f s sf

f

( )

( )( ) ( () )( ) 2 )

( )

2

1

3 3

2

For the external portion i g of the error function, we find similarly: 
27

(A.3.7) 

( )( ) .

1( ) ( )

2

1

( )

2

1 2 2

















+ ∇ − + ∇ 















= ∇ −

















= ∇ ∇ = ∇ − + −

i

T i

i i

i

i

i

T

i

T

i

i

i

i

i

T

XYZ i

T

i i

q

Rp d

q

Rp

q

R

d

q

Rp

q

R

g g

i i i i i i i i i

i i i i i i

ξ η ς ξ η ς ξ η ς

ξ η ς ξ η ς

U Λ Λ

Γ U Λ

 . 

By (A.3.4), 

T

i

i

i

T

q

Rp

q

R

i i i

= − 3 U 















∇ξ η ς − , 

( )

T

i

T

i

T i i

i

i

i

i

T

i i

T

i i

T i

q

R X Y Z

q

R q p

q

Rp

q

R

q

p Rp

q

R

q

Rp

i i i i i i i i i

U U U













 + +

= 













 +

= 















= +

= ∇ + ∇ 















∇

3

2 2 2

3

2 2

3

2

( ) ( )

1

ξ η ς ξ η ς ξ η ς

so that – in matrix notation--, 

 .

( )

( ) ( )

2

1

3

2 2 2

3

2 T

i

i i

i

i

i

T T

i

i

i i

q

R X Y Z

d

q

Rp

q

Rp

g

i i i

Γ U U Λ Λ U













 + +

+















+ ∇ − 















= − ξ η ς

and finally, 

(A.3.8) Γ UU ΛU I

















+ − 













 + +

+ 















= − i

i

T i

i

i

T

i

i

i i d

q

Rp

q

R X Y Z

q

Rp

g 3

2 2 2

3

2

( )

( )

2

1

 . 

A.4 Gradients and Hessians of the Orthogonal Error Function 

We repeat the definitions of Chapter 4. The error function for orthogonal fitting of a sphere is 

given as 

 

2 2

1

2

1

2 E h (w R) wi

2Rwi R

n

i

i

n

i

orth = ∑ i = ∑ − = − +

= =

where 
28

 hi = wi − R 

and 

 

2 2 2 2 2 2 2 2

i

( i) ( i) ( i) (X Y Z ) 2di pi di w = X − x + Y − y + Z − z = + + − + 

 with the notations (3.1.1) 

 

 i i i i i ix = d cos , φ cosθ = d ξ i i i i i i

y = d sin φ cos θ = d η , i i i i i

z = d sin , θ = d ξ

and the definition (3.1.2) of the auxiliary quantity i p . The vector 

. 

 i i

i

i

i

i

i

i

i

i d d

Z

Y

X

Z z

Y y

X x

W = U − Λ





























−





























=





























−

−

−

=

ς

η

ξ

 

will play the key role. Indeed, all the following gradients are multiples of Wi

: 

(A.4.1) XYZwi 2Wi

2 ∇ =

 i

i

XYZ i

i

XYZ i XYZ i

w

w

w

w w W

1 1

2

2 1 2 ∇ = ∇ = ∇ = , 

 i

i

XYZ i

i i

XYZ

w

w

w w

2 3 W

1 1 1

∇ = − ∇ = − . 

 

Thus 

(A.4.2) i

i

XYZ i XYZ i

w

h w W

1

∇ = ∇ = , 

and in view of 

 i

i

i

XYZ i XYZ i i XYZ i

w

w R

h w R w R w W

−

∇ = ∇ ( − ) = (2 − )∇ = 2

2 2

 , 

it follows that 
29

(A.4.3) i

i

XYZ i

w

R

h



W













∇ = 1−

2

1 2

 . 

 

Te Hessians of the quantities considered below are linear combinations of the two matrices 

 I = identity, and [ ] T

i i i i i

i

i

i

X x Y y Z z

Z z

Y y

X x

− − − = W W





























−

−

−

. 

In particular, 

 

(A.4.4) H H ( ) 2 H H H ( ) 2I

2 2 2 2 2 2 2 2

XYZwi = XYZ X + Y + Z − di XYZ pi + XYZ di = XYZ X + Y + Z =

and by (A.1.2), 

 

.

1 1

2

1

2

1

2 2

1

4

1 1

2

1 1

4

1

3

3

2 2 2

3

W W I

H H W W I

i

T

i i

i

i

T

i i

i

XYZ i

i

i

T

XYZ i XYZ

i

XYZ i

w w

w w

w

w

w w

w

w

= − +

= − ∇ ∇ + = − +

 

As hi = wi − R and 2 2 2

hi = wi − 2Rwi + R , we have 

 H H WW I

i

T

i i

i

XYZ i XYZ i

w w

h w

1 1

3

= = − + . 

and 

 H H H I WW I

i

T

i i

i

XYZ i XYZ i XYZ i

w

R

w

R

h w R w

2 2

2 2 3

2 2

= − = + − . 

so that 

(A.4.5) T

i i

i i

XYZ i

w

R

w

R

H h I

3 WW

2

1

2

1

+ 















= − . 

A.5 Mixed Derivatives of the Orthogonal Error Function 
30

First, we differentiate with respect to the range variable i d , starting with the key quantity wi : 

(A.5.1) ( 2 ) (2 )

2 2 2 2

2

i i i i i

i i

i X Y Z d p d p d

d d

w

+ + − + = − −

∂

∂

=

∂

∂

 ( )

1

( )(2 )

1

2

1 1

2

1

2

i i

i

i i

i i

i

i i

i

p d

w

p d

d w

w

d w

w

= − − = − −

∂

∂

=

∂

∂

 ( )

1

( )(1 )

1 1 1

2 3 3 i i

i

i i

i i i i i

p d

w

p d

d w

w

d w w

= − − − = −

∂

∂

= −

∂

∂

 . 

Thus 

 ( ) (2 ) 2 ( 2 1 ( )

2

2

i i

i

i i

i

i

i

i

i i

i i

i

p d

w

R

p d

w

w R

d

w

w R w R

d d

h

− 















− = − −

−

= −

∂

∂

− = −

∂

∂

=

∂

∂

 

and 

 XYZ i

i

i i

i

XYZ

i

i

XYZ p

w

R

p d

w

R

d

h

∇















− − − 















= − − ∇

∂

∂

∇ ( ) 1

1

2

1

2

so that, finally, by (A.4.1) and (3.2.1), 

(A.5.2) i

i

i i i

i i

i

XYZ

w

R

p d

w

R

d

h W Λ















= − − − −

∂

∂

∇ ( ) 1

2

1

3

2

. 

We move to differentiation with respect to the bearing variables φi θi

, . Again we aim to apply 

the differential matrix operator 

 [ ] XYZ

T

i

i i i

Γ = ∇ξ η ς ∇ , 

here to the individual orthogonal error 2

i h -- expressed in matrix notation --: 

















+ ∇ − 















= ∇ − 















 = ∇ −











= ∇ ∇

i

T

i

T

i

T

XYZ i

T

i i

w

R

w

R

w

R

h h

i i i i i i i i i i i i

1 ( ) 1 1

2

1

( )

2

1 2 2 Γ ξ η ς ξ η ς W ξ η ς W W ξ η ς

. 

By (A.3.3), 
31

 W Λ I i i

T

i i

T

d d

i i i i i i

∇ξ η ς = − ∇ξ η ς = − , 

and by (A.3.4), 

 ( )

T

i i

T

i i i i

T T

w X Y Z d p d d p d

i i i i i i i i i

 ( ) 2 2 2 U

2 2 2 2 2 ∇ξ η ς i = ∇ξ η ς + + − + = − ∇ξ η ς = −

 

T

i

T T T i

w

d

w

w

w w i i i i i i i i i

∇ = ∇ = ∇ = − U

2

i

2

i

2

i i

1

2

1

 ξ η ς ξ η ς ξ η ς

 

T

i

T i

i

T T i

w

d

w

d

w w w

i i i i i i

U 3 U

i

i

i i

1

w

1 1

 =















∇ξ η ς = − ∇ξ η ς = − −

so that by (A.4.1), 

 

T

i

i

i

i i i

w

Rd

w

R

Γ h d I

3 WU 2

( ) 1

2

1

+ 















= − − . 

Taking into account the orthogonality relations (3.2.10), this matrix can again be symmetrized by 

substituting T

i

T T W = U + Λ for T U in the above expression, yielding the matrix 

 

T

i

i

i

i i

S

i

w

Rd

w

R

Γ h d I

3 WW 2

( ) 1

2

1

+ 















= − − . 

Acknowledgement 

 

We are grateful to our NIST colleagues James Filliben, Alan Heckert, and Joan Rosenblatt for 

discussions and material about the subject of variance propagation. 

REFERENCES 

[1] W. C. Stone, M. Juberts, N. Dagalakis, J. Stone, and J. Gorman, "Performance Analysis 

of Next-Generation LADAR for Manufacturing, Construction, and Mobility," NISTIR 

vol. 7117, 2004. 

[2] D. E. Gilsinn, C. Witzgall, G. S. Cheok, and A. Lytle, "Construction Object Identification 

from LADAR Scans: An Experimental Study Using I-Beams," NISTIR, vol. 7286, 2005. 
32

[3] P. J. Besl and N. D. McKay, "A method for registration of 3-D shapes," IEEE Trans. 

Pattern Analysis and Machine Intelligence, vol. 14, pp. 239-256, 1992. 

[4] R. O. Duda and P. E. Hart, "Use of Hough transformation to detect lines and curve in 

pictures," Comm. ACM, vol. 15, pp. 11 - 15, 1972. 

[5] K. Kanatani, Statistical Optimization for Geometric Computation. New York: Elsevier, 

1996. 

[6] A. C. Davidson, Statistical Models. Cambridge: Cambridge University Press, 2003. 

[7] D. A. Freedman, Statistical Models, Theory and Practice, revised ed. Cambridge: 

Cambridge University Press, 2009. 

[8] P. T. Boggs, J. R. Donaldson, R. B. Schnabel, and C. H. Spiegelman, "A computational 

examination of orthogonal distance regression," Journal of Econometrics, vol. 38, pp. 

169 - 201, 1988. 

[9] P. T. Boggs, J. R. Donaldson, R. H. Byrd, and R. B. Schnabel, "ALGORITHM 676 

ODRPACK: Software for Weighted Orthogonal Distance Regression," ACM 

Transactions on Mathematical Software, vol. 15, pp. 348 - 364, 1989. 

[10] W. Gander, G. H. Golub, and R. Strebel, "Least-squares fitting of circles and ellipses," 

BIT, vol. 34, pp. 558-578, 1994. 

[11] C. M. Shakarji, "Least-Squares Fitting Algorithms of the NIST Algorithm Testing 

System," Journal of Research of NIST, vol. 103, pp. 633-640, 1998. 

[12] S. J. Ahn, W. Rauh, and H. J. Warnecke, "Least-squares orthogonal distances fitting of 

circle, sphere, ellipse, hyperbola, and parabola," Pattern Recognition, vol. 34, pp. 2283-

2303, 2001. 

[13] I. D. Coope, "Circle Fitting by Linear and Nonlinear Least Squares," J. Optimization 

Theory Applicat., vol. 76, pp. 381-388, 1993. 

[14] B. A. Jones and R. B. Schnabel, "A comparison of two sphere fitting methods.," in IEEE 

Proceedings of the Instrumentation and Measurement Technology, 1986. 

[15] I. Kasa, "A Circle Fitting Procedure and Its Error Analysis," IEEE Trans. Instrum. Meas., 

vol. 25, pp. 8-14, 1976. 

[16] Y. Nievergelt, "Computing circles and spheres of arithmetic least squares," Comput. 

Phys. Commun., vol. 81, pp. 343-350, 1994. 

[17] M. Renault, "Fitting Circles and Ellipses to Data Using the Least-Squares Method," in 

http://www.math.temple.edu/~renault/ellipses.html. 

[18] H. Spath, "Least-Square Fitting with Spheres," J. Optimization Theory Applicat., vol. 96, 

pp. 191-199, 1998. 

[19] D. Umbach and K. N. Jones, "A Few Methods for Fitting Circles to Data," IEEE Trans. 

Instrum. Meas. , vol. 52, pp. 1881-1885, 2003. 

[20] C. Witzgall, G. S. Cheok, and A. J. Kearsley, "Recovering Circles and Spheres from 

Point Data," in Perspectives In Operations Research, F. B. Alt, M. C. Fu, and B. L. 

Golden, Eds. New York: Springer, 2006, pp. 393-413.

[21] M. Franaszek, G. S. Cheok, K. S. Saidi, and C. Witzgall, "Fitting Spheres to Range Data 

from 3D Imaging Systems," IEEE Trans. Instrum. Meas. , 2009 - in print. 

[22] M. Franaszek, G. S. Cheok, K. S. Saidi, and C. Witzgall, "Sphere Fitting as a Measure of 

3D Imaging Systems Performance," in preparation, 2009. 

[23] C. Witzgall and M. Franaszek, "Discussion of Variance Calculation for Linear and 

Nonlinear Regression," in preparation, 2009. 
33

[24] GUM, "Evaluation of measurement data - Guide to the expression of uncertainty in 

measurement (http://www.bipm.org/en/publications/guides/gum.html)," JCGM 2008. 

[25] H. Ku, "Notes on the Use of Propagation of Error Formulas," Journal of Research 

National Bureau of Standards, vol. 70C, pp. 331 - 341, 1966. 

[26] D. M. Himmelblau, Process Analysis by Statistical Methods. New York: Wiley, 1970. 

[27] N. R. Draper and H. Smith, Applied Regression Analysis, 2nd ed. New York: Wiley, 

1981. 

[28] C. G. Broyden, "The Convergence of a Class of Double-rank Minimization Algorithms: 

2. The New Algorithm " IMA J. Appl. Math., vol. 6, pp. 222 - 231, 1970.