Document [original]

mathematics

Article

Weighted Total Least Squares (WTLS) Solutions for

Straight Line Fitting to 3D Point Data

Georgios Malissiovas 1,* , Frank Neitzel 1, Sven Weisbrich 1and Svetozar Petrovic 1,2

1Institute of Geodesy and Geoinformation Science, Technische Universität Berlin, Strasse des 17. Juni 135,

10623 Berlin, Germany; [email protected] (F.N.); [email protected] (S.W.);

svetozar[email protected] (S.P.)

2GFZ German Research Centre for Geosciences, Section 1.2: Global Geomonitoring and Gravity Field,

Telegrafenberg, 14473 Potsdam, Germany

*Correspondence: geor[email protected]

Received: 14 July 2020; Accepted: 25 August 2020; Published: 29 August 2020





Abstract:

In this contribution the fitting of a straight line to 3D point data is considered, with Cartesian

coordinates

as observations subject to random errors. A direct solution for the case of equally

weighted and uncorrelated coordinate components was already presented almost forty years ago.

For more general weighting cases, iterative algorithms, e.g., by means of an iteratively linearized

Gauss–Helmert (GH) model, have been proposed in the literature. In this investigation, a new direct

solution for the case of pointwise weights is derived. In the terminology of total least squares (TLS),

this solution is a direct weighted total least squares (WTLS) approach. For the most general weighting

case, considering a full dispersion matrix of the observations that can even be singular to some extent,

a new iterative solution based on the ordinary iteration method is developed. The latter is a new

iterative WTLS algorithm, since no linearization of the problem by Taylor series is performed at any

step. Using a numerical example it is demonstrated how the newly developed WTLS approaches

can be applied for 3D straight line fitting considering different weighting cases. The solutions are

compared with results from the literature and with those obtained from an iteratively linearized

GH model.

Keywords:

3D straight line fitting; total least squares (TLS); weighted total least squares (WTLS);

nonlinear least squares adjustment; direct solution; singular dispersion matrix; laser scanning data

1. Introduction

Modern geodetic instruments, such as terrestrial laser scanners, provide the user directly with

3D coordinates in a Cartesian coordinate system. However, in most cases these 3D point data are not

the final result. For an analysis of the recorded data or for a representation using computer-aided design

(CAD), a line, curve or surface approximation with a continuous mathematical function is required.

In this contribution the fitting of a spatial straight line is discussed considering the coordinate

components

of each point

as observations subject to random errors, which results in

a nonlinear adjustment problem. An elegant direct least squares solution for the case of equally

weighted and uncorrelated observations has already been presented in 1982 by Joviˇci´c et al. [

Unfortunately, this article was very rarely cited in subsequent publications, which is probably due

to the fact that it was written in Croatian language. Similar least squares solutions, direct as well,

have been published by Kahn [

] and Drixler ([

], pp. 46–47) some years later. In these contributions,

it was shown that the problem of fitting a straight line to 3D point data can be transformed into an

eigenvalue problem. An iterative least squares solution for fitting a straight line to equally weighted

Mathematics 2020,8, 1450; doi:10.3390/math8091450 www.mdpi.com/journal/mathematics

Mathematics 2020,8, 1450 2 of 19

and uncorrelated 3D points has been presented by Späth [

], by minimizing the sum of squared

orthogonal distances of the observed points to the requested straight line.

For more general weighting schemes iterative least squares solutions have been presented

by Kupferer [

] and Snow and Schaffrin [

]. In both contributions various nonlinear functional

models were introduced and tested and the Gauss-Newton approach has been employed for

an iterative least squares solution, which involves the linearization of the functional model

to solve the adjustment problem within the Gauss–Markov (GM) or the Gauss–Helmert (GH)

models. The linearization of originally nonlinear functional models is a very popular approach

in adjustment calculation,

where the solution

is then determined by iteratively linearized GM or GH

models,

see e.g., the textbooks

by Ghilani [

], Niemeier [

] or Perovic [

]. Pitfalls to be avoided in

the iterative adjustment of nonlinear problems have been pointed out by Pope [

] and Lenzmann

and Lenzmann [11].

Fitting a straight line to 3D point data can also be considered as an adjustment problem of

type total least squares (TLS) for an errors-in-variables (EIV) model, as already pointed out by

Snow and Schaffrin [

]. For some weighting cases, problems expressed within the EIV model can

have a direct solution using singular value decomposition (SVD). This solution is known under

the name Total Least Squares (TLS) and has been firstly proposed by Golub and Van Loan [

Van Huffel and Vandewalle [13]

or Van Huffel and Vandewalle ([

], p. 33 ff.). The TLS solution is

obtained by computing the roots of a polynomial, i.e., by solving the characteristic equation of the

eigenvalues, which is identical in the case of fitting a straight line to 3D point data to the solutions

of Joviˇci´c et al. [

], Kahn [

] and Drixler ([

], pp. 46–47). This is something that has been observed

by Malissiovas et al. [

], where a relationship has been presented between TLS and the direct

least squares solution of the same adjustment problem, while postulating uncorrelated and equally

weighted observations.

To involve more general weight matrices in the adjustment procedure, iterative algorithms have

been presented in the TLS literature without linearizing the underlying problem by Taylor series at any

step of the solution process. These are algorithmic approaches known as weighted total least squares

(WTLS), presented e.g., by Schaffrin and Wieser [

], Shen et al. [

] or Amiri and Jazaeri [

]. A good

overview of such algorithms, as well as alternative solution strategies can be found in the dissertation

theses of Snow [

] and Malissiovas [

]. An attempt to find a WTLS solution for straight line fitting to

3D point data was made by Guo et al. [21].

To avoid confusion, it is to clarify that the terms TLS and WTLS refer to algorithmic approaches

for obtaining a least squares solution, which is either direct or iterative but without linearizing the

problem by Taylor series at any step. This follows the statement of Snow ([

], p. 7), that “the terms

total least-squares (TLS) and TLS solution [. ..] will mean the least-squares solution within the EIV

model without linearization”. Of course, a solution within the GH model is more general in the sense

that it can be utilized to solve any nonlinear adjustment problem, while TLS and WTLS algorithms can

treat only a certain class of nonlinear adjustment problems. This has been firstly discussed by Neitzel

and Petrovic [

] and Neitzel [

], who showed that the TLS estimate within an EIV model can be

identified as a special case of the method of least squares within an iteratively linearized GH model.

To the extent of our knowledge, a WTLS algorithm for fitting a straight line to 3D point data has

not been presented yet. Therefore, in this study we derive two novel WTLS algorithms for the discussed

adjustment problem considering two different cases of stochastic models:

(i)

pointwise weights, i.e., coordinate components with same precision for each point and no

correlations between them,

(ii) general weights, i.e., correlated coordinate components of individual precision including singular

dispersion matrices.

The adjustment problem resulting from case (i) can still be solved directly, i.e., a direct WTLS

solution, presented in Section 2.1 of this paper. For case (ii) an iterative solution without linearizing

Mathematics 2020,8, 1450 3 of 19

the problem by Taylor series is derived in Section 2.2, i.e., an iterative WTLS solution. Both solutions

are based on the work of Malissiovas [

], where similar algorithms have been presented for the

solution of other typical geodetic tasks, such as straight line fitting to 2D point data, plane fitting to 3D

point data and 2D similarity coordinate transformation. The WTLS solution for straight line fitting

to 3D point data will be derived from a geodetic point of view by means of introducing residuals for

all observations, formulating appropriate condition and constraint equations, setting up and solving

the resulting normal equation systems.

2. Straight Line Fitting to 3D Point Data

A straight line in 3D space can be expressed by

y−y0

a=x−x0

b=z−z0

c, (1)

as explained e.g., in the handbook of mathematics by Bronhstein et al. ([

], p. 217). This equation

defines a line that passes through a point with coordinates

and

and is parallel to a direction

vector with components

and

. Since the number of unknown parameters is six, two additional

constraints between the unknown parameters have to be taken into account for a solution, as Snow

and Schaffrin [

] have clearly explained. Proper constraints will be selected at a later point of

the derivations.

Considering an overdetermined configuration for which we want to obtain a least squares solution,

we introduce for all points

residuals

vxi

vyi

vzi

for the observations

assuming that they are

normally distributed with

v∼(0,σ2

0QLL)

. Here the 1

vector

contains the residuals

vxi

vyi

vzi

QLL

is the 3

n×

corresponding cofactor matrix and

σ2

the theoretical variance factor. Based on (1),

the nonlinear conditions equations

a(xi+vxi−x0)−b(yi+vyi−y0) = 0,

b(zi+vzi−z0)−c(xi+vxi−x0) = 0,

c(yi+vyi−y0)−a(zi+vzi−z0) = 0,

(2)

can be formulated for each point

, with

. . .

and

being the number of observed points.

A functional model for this problem can be expressed by two of these three nonlinear condition

equations per observed point. Using all three condition equations for solving the problem would lead

to a singular normal matrix due to linearly dependent normal equations.

2.1. Direct Total Least Squares Solution for Equally Weighted and Uncorrelated Observations

Considering all coordinate components

for all points

as equally weighted and

uncorrelated observations with

pyi=pxi=pzi=1∀i, (3)

a least squares solution can be derived by minimizing the objective function

Ω=

∑

i=1

xi+v2

yi+v2

zi→min. (4)

A direct least squares solution of this problem, respectively a TLS solution, has been presented by

Malissiovas et al. [

] and Malissiovas [

]. According to these investigations, equally weighted

residuals correspond to the normal distances

i=v2

xi+v2

yi+v2

zi, (5)

Mathematics 2020,8, 1450 4 of 19

as deviation measures between the observed 3D point cloud and the straight line to be computed.

Thus, the objective function can be written equivalently as

Ω=

∑

i=1

xi+v2

yi+v2

zi=

∑

i=1

i→min. (6)

An expression for the squared normal distances

D2=[a(x−x0)−b(y−y0)]2+[b(z−z0)−c(x−x0)]2+[c(y−y0)−a(z−z0)]2

a2+b2+c2,(7)

can be found in the handbook of Bronhstein et al. ([24], p. 218).

An appropriate parameterization of the problem involves the substitution of the unknown

parameters

and

with the coordinates of the center of mass of the observed 3D point data.

A proof for this parameter replacement has been given by Joviˇci´c et al. [

] and concerns only the case of

equally weighted and uncorrelated observations. A solution for the line parameters can be computed

by minimizing the objective function (6), under the constraint

a2+b2+c2=1, (8)

or by searching for stationary points of the Lagrange function

∑

i=1

i−k(a2+b2+c2−1), (9)

with

denoting the Lagrange multiplier. The line parameters can be computed directly from the normal

equations, either by solving a characteristic cubic equation or an eigenvalue problem. A detailed

explanation of this approach was given by Malissiovas ([

], p. 74 ff.). In the following subsections

we will give up the restriction of equally weighted and uncorrelated coordinate components and derive

least squares solutions considering more general weighting schemes for the observations.

2.2. Direct Weighted Total Least Squares Solution

In this case we consider the coordinate components

of each point

to be uncorrelated

and of equal precision with

σyi=σxi=σzi=σi∀i, (10)

yielding the corresponding (pointwise) weights

pi=1

σ2

∀i. (11)

A least squares solution of the problem can be found by minimizing the objective function

Ω=

∑

i=1

pxiv2

xi+pyiv2

yi+pziv2

zi→min

∑

i=1

piv2

xi+v2

yi+v2

zi→min.

(12)

Taking into account the relation of the squared residuals with the normal distances of Equation (5), it is

possible to express the objective function as

Ω=

∑

i=1

piD2

i→min. (13)

Mathematics 2020,8, 1450 5 of 19

Using the constraint (8), the expression of the normal distances can be written as

i=[a(xi−x0)−b(yi−y0)]2+[b(zi−z0)−c(xi−x0)]2+[c(yi−y0)−a(zi−z0)]2.(14)

A further simplification of the problem is possible, by replacing the unknown parameters

and

with the coordinates of the weighted center of mass of the 3D point data

y0=

∑

i=1

piyi

∑

i=1

,x0=

∑

i=1

pixi

∑

i=1

,z0=

∑

i=1

pizi

∑

i=1

. (15)

It can be proven that this parameterization is allowed and possible, following the same line of thinking

as in the study of Joviˇci´c et al. [

]. Therefore, for this weighted case we can write the normal distances as

i=ax0

i−by0

i2+bz0

i−cx0

i2+cy0

i−az0

i2,(16)

with the coordinates reduced to the weighted center of mass of the observed 3D point data

i=yi−

∑

i=1

piyi

∑

i=1

,x0

i=xi−

∑

i=1

pixi

∑

i=1

,z0

i=zi−

∑

i=1

pizi

∑

i=1

. (17)

Searching for stationary points of the Lagrange function

∑

i=1

i−k(a2+b2+c2−1)

∑

i=1ax0

i−by0

i2+bz0

i−cx0

i2+cy0

i−az0

i2−k(a2+b2+c2−1)

(18)

leads to the system of normal equations

∂K

∂a=2a n

∑

i=1

pix0

∑

i=1

piz0

2−k!−2b

∑

i=1

piy0

ix0

i−2c

∑

i=1

piy0

iz0

i=0, (19)

∂K

∂b=2b n

∑

i=1

piy0

∑

i=1

piz0

2−k!−2a

∑

i=1

piy0

ix0

i−2c

∑

i=1

pix0

iz0

i=0, (20)

∂K

∂c=2c n

∑

i=1

piy0

∑

i=1

pix0

2−k!−2a

∑

i=1

piy0

iz0

i−2b

∑

i=1

pix0

iz0

i=0, (21)

and ∂K

∂k=−a2+b2+c2−1=0. (22)

Equations (19) to (21) are a homogeneous system of equations which are linear in the unknown

line parameters a,band c, and has a nontrivial solution if



(d1−k)d4d5

d4(d2−k)d6

d5d6(d3−k)

=0, (23)

Mathematics 2020,8, 1450 6 of 19

with the elements of the determinant

d1=

∑

i=1

pix0

∑

i=1

piz0

2,d2=

∑

i=1

piy0

∑

i=1

piz0

2,d3=

∑

i=1

piy0

∑

i=1

pix0

d4=−

∑

i=1

piy0

ix0

i,d5=−

∑

i=1

piy0

iz0

i,d6=−

∑

i=1

pix0

iz0

(24)

Equation (23) is a cubic characteristic equation with the unknown parameter

and a solution is

given by Bronhstein et al. ([

], p. 63). The unknown line parameters

and

can be computed either

by substituting

kmin

into Equations (19)–(21) under the constraint (22) or by transforming the equation

system and solving an eigenvalue problem, i.e., the direct WTLS solution.

2.3. Iterative Weighted Total Least Squares Solution

In this case, we consider the fact that the coordinate components

of each point

are not

original observations. They are either

(i)

derived from single point determinations, e.g., using polar elementary measurements of slope

distance

ρi

, direction

φi

and tilt angle

θi

, e.g., from measurements with a terrestrial laser scanner,

so that the coordinates result from the functional relationship

xi=ρicos φisin θi,

yi=ρisin φicos θi,

zi=ρicos θi.

(25)

Using the standard deviations

σρi

σφi

and

σθi

of the elementary measurements, a 3

variance-covariance matrix can be provided for the coordinate components

of each

point

by variance-covariance propagation based on (25). Covariances among the

different

points do not occur in this case.

(ii)

derived from a least squares adjustment of an overdetermined geodetic network.

From the solution within the GM or GH model, the 3

n×

variance-covariance matrix of

the coordinate components

of all points

can be obtained. This matrix is a full matrix,

considering also covariances among the

different points. It may even have a rank deficiency in

case the coordinates were determined by a free network adjustment.

The variances and covariances from (i) or (ii) can then be assembled in a dispersion matrix

QLL =





Qxx Qxy Qxz

Qyx Qyy Qyz

Qzx Qzy Qzz





, (26)

for the coordinate components as “derived observations”. In case of a non-singular dispersion matrix,

it is possible to compute the respective weight matrix

P=





Pxx Pxy Pxz

Pyx Pyy Pyz

Pzx Pzy Pzz





=Q−1

LL . (27)

Starting with the functional model (2), we choose to work with the nonlinear condition equations

c(yi+vyi−y0)−a(zi+vzi−z0) = 0,

c(xi+vxi−x0)−b(zi+vzi−z0) = 0.

(28)

Mathematics 2020,8, 1450 7 of 19

Two additional constraints must be taken into account in this case. For example,

Snow and Schaffrin [6]

proposed the constraints

a2+b2+c2=1,

ay0+bx0+cz0=0.

(29)

A further development of an algorithmic solution using these constraints is possible.

However, we choose

to parameterize the functional model so that an additional constraint is

unnecessary and the equations become simpler. Therefore, we choose to fix one coordinate of the point

on the line

z0=

∑

i=1

n=z, (30)

as well as one of the unknown components of the direction vector

c=1. (31)

This simplification of the functional model has been used by Boroviˇcka et al. [

] for a similar practical

example, this of estimating the trajectory of a meteor. As already mentioned in that work, the resulting

direction vector of the straight line can be transformed to a unit vector afterwards, i.e., a solution using

the constraint a2+b2+c2=1. The simplified functional model can be written in vector notation as

(yc+vy−ey0)−a(zc+vz−ez) = 0,

(xc+vx−ex0)−b(zc+vz−ez) = 0,

(32)

with the vectors

and

containing the coordinates of the observed points and vectors

and

the respective residuals. Vector

is a vector of ones, with length equal to the number of

observed points n.

A weighted least squares solution of this problem can be derived by minimizing the

objective function

Ω=vT

xPxxvx+vT

yPyyvy+vT

zPzzvz+2vT

xPxyvy+2vT

xPxzvz+2vT

yPyzvz(33)

or by searching for stationary points of the Lagrange function

K=Ω−2kT

1(yc+vy−ey0)−a(zc+vz−ez)

−2kT

2[(xc+vx−ex0)−b(zc+vz−ez)](34)

with two distinct vectors

and

of Lagrange multipliers. Taking the partial derivatives with respect

to all unknowns and setting them to zero results in the system of normal equations

∂K

∂vT

=Pxxvx+Pxyvy+Pxzvz−k2=0,(35)

∂K

∂vT

=Pyyvy+Pyxvx+Pyzvz−k1=0,(36)

∂K

∂vT

=Pyyvy+Pxzvx+Pyzvy+ak1+bk2=0,(37)

Mathematics 2020,8, 1450 8 of 19

∂K

∂kT

= (yc+vy−ey0)−a(zc+vz−ez) = 0, (38)

∂K

∂kT

= (xc+vx−ex0)−a(zc+vz−ez) = 0, (39)

∂K

∂a=kT

1(zc+vz−ez)=0, (40)

∂K

∂b=kT

2(zc+vz−ez)=0, (41)

∂K

∂y0

=eTk1=0, (42)

∂K

∂x0

=eTk2=0. (43)

To derive explicit expressions for the residual vectors, we write Equations (35)–(37) as







Pxx Pxy Pxz

Pyx Pyy Pyz

Pzx Pzy Pzz















=





−ak1−bk2





, (44)

which leads to the solution for the residual vectors











=





Pxx Pxy Pxz

Pyx Pyy Pyz

Pzx Pzy Pzz







−1





−ak1−bk2







=





Qxx Qxy Qxz

Qyx Qyy Qyz

Qzx Qzy Qzz











−ak1−bk2





,

(45)

or explicitly

vx=Qxy −aQxzk1+(Qxx −bQxz)k2, (46)

vy=Qyy −aQyzk1+Qxy −bQyzk2, (47)

vz=Qyz −aQzk1+(Qxz −bQzz)k2. (48)

Inserting these expressions for the residual vectors into Equations (38) and (39) yields

W1k1+W2k2+yc−ey0−azc=0, (49)

W2k1+W3k2+xc−ex0−bzc=0, (50)

Mathematics 2020,8, 1450 9 of 19

with the auxiliary matrices W1,W2and W3as

W1=Qy−2aQyz +a2Qz, (51)

W2=Qxy −aQxz −bQyz +abQz, (52)

W3=Qx−2bQxz +b2Qz, (53)

Using Equations (49), (50) and (40)–(43), the nonlinear equation system

W1k1+W2k2−ey0−azc=−yc,

W2k1+W3k2−ex0−bzc=−xc,

−(zc+vz−ez)Tk1=0,

−(zc+vz−ez)Tk2=0,

−eTk1=0,

−eTk2=0

(54)

can be set up. To express this system of equations into a block matrix form and to obtain a symmetrical

matrix of normal equations in the following, it is advantageous to add the terms

−a(vz−ez)

and −b(vz−ez)to both sides of the first two equations, leading to the system of equations

W1k1+W2k2−ey0−a(zc+vz−ez)=−yc−a(vz−ez),

W2k1+W3k2−ex0−b(zc+vz−ez)=−xc−b(vz−ez),

−(zc+vz−ez)Tk1=0,

−(zc+vz−ez)Tk2=0,

−eTk1=0,

−eTk2=0.

(55)

Arranging the unknowns in a vector













, (56)

the equation system (55) can be written as

NX =n. (57)

Mathematics 2020,8, 1450 10 of 19

with the matrix of normal equations

N="W A

AT0#, (58)

constructed using the matrices

W="W1W2

W2W3#(59)

and

AT=





−(zc+vz−ez)0

0−(zc+vz−ez)

−e 0

0−e







. (60)

The right-hand side of the equation system (57) reads







−yc−a(zc+vz−ez)

−xc−b(zc+vz−ez)







(61)

It is important to point out that matrix

and vector

contain the unknown parameters

and

in their entries. Therefore, to express and solve these normal equations that live in the “nonlinear

world” with the help of vectors and matrices (that only exist in the “linear world”), appropriate

approximate values

and

have to be introduced for all auxiliary matrices required for the setup

of matrix Nand vector n. The solution for the vector of unknowns can be computed by

X=N−1n, (62)

without the need of a linearization by Taylor series at any step of the calculation process. The WTLS

solution for the line parameters can be computed following the ordinary iteration method, as it is

explained for example by Bronhstein ([

], p. 896). Therefore, the solutions

and

˜vz

, after stripping

them of their random character, are to be substituted as new approximate values as long as necessary

until a sensibly chosen break-off condition is met. For example, the maximum absolute difference

between two consecutive solutions must be smaller than a predefined threshold

, which in this

problem can be formulated as

max "a0

b0#−"ˆ

b#

≤e, (63)

with

|·|

denoting the absolute value. The predicted residual vectors

˜vx

˜vy

˜vz

can be computed

from Equations (46)–(48). The iterative procedure for the presented WTLS solution can be found in

Algorithm 1.

Mathematics 2020,8, 1450 11 of 19

Algorithm 1 Iterative WTLS solution

Choose approximate values for a0,b0and v0

Define parameter c=1.

Set threshold efor the break-off condition of the iteration process.

Set parameter da=db=∞, for entering the iteration process.

while da>eand db>edo

Compute matrices Wand A.

Build matrix Nand vector n.

Estimate the vector of unknowns ˆ

Compute the residual vector ˜vz.

Compute parameters da=|ˆ

a−a0|and db=|ˆ

b−b0|.

Update the approximate values with the estimated ones, with a0=ˆ

a,b0=ˆ

band v0

z=˜vz.

end while

return ˆ

aand ˆ

b, with c=1.

After computing the line parameters

and

and putting them into a vector, we can easily scale

it into a unit vector by dividing each component with the length of the vector











=





























(64)

with

||·||

denoting the Euclidean norm. The derived parameters

and

refer to the normalized

components of a unit vector, that is parallel to the requested line, with

n+b2

n+c2

n=1 (65)

2.4. WTLS Solution with Singular Dispersion Matrices

The algorithmic approach presented in Section 2.3 can also cover cases when dispersion matrices

are singular. Such a solution depends on the inversion of matrix

(58), which depends on the rank

deficiency of matrix

(59). Following the argumentation of Malissiovas ([

], p. 112), a criterion that

ensures a unique solution of the problem can be described in this case by

rank ([W|A]) =2n, (66)

with

rank of

W≤

, with 2

= number of condition equations, since for each of the observed

points

two condition equations from Equation (2) are taken into account;

- rank of A=m, with m= number of unknown parameters.

In cases of singular dispersion matrices, the rank of matrix

will be smaller than 2

. A unique

solution will exist if the rank of the augmented matrix

[W|A]

is equal to the number of condition

equations 2

of the problem. It is important to mention that the developed idea is based on

the Neitzel–Schaffrin criterion, which has been firstly proposed by Neitzel and Schaffrin [

particularly for a solution of an adjustment problem within the GH model when singular dispersion

matrices must be employed.

Mathematics 2020,8, 1450 12 of 19

2.5. A Posteriori Error Estimation

In this section we want to determine the variance-covariance matrix of the estimated parameters.

The following derivations can be used for computing the a posteriori stochastic information for

all weighted cases discussed in this investigation, i.e., the direct and the iterative WTLS solutions.

Therefore, we will employ the fundamental idea of variance-covariance propagation. This is a

standard procedure explained by many authors, like for example in the textbooks of Wells and

Krakiwsky ([

], p. 20), Mikhail ([

], p. 76 ff.) or Niemeier ([

], p. 51 ff.) and has been further employed

in the GM and the GH model, so that the a posteriori stochastic results can be computed using directly

the developed matrices from each model. A detailed explanation is given by Malissiovas ([

], p. 31 ff.).

As we have already mentioned in the introduction of this article, a TLS, respectively a WTLS

solution, can be regarded as a special case of a least squares solution within the GH model.

From the presented

WTLS algorithm we observe that the derived matrix of normal Equation (58)

is equal to the matrix if it was computed within the GH model. Therefore, it is possible to compute

N−1="Q11 Q12

Q21 Q22 #(67)

and extract the dispersion matrix for the unknown parameters from

Qˆ

xˆ

x=−Q22. (68)

The a posteriori variance factor is

σ2

0=vTPv

r(69)

with vector

holding all residuals and

denoting the redundancy of the adjustment problem. In case

of a singular dispersion matrix, it is not possible to compute the weight matrix

, as in Equation (27).

Therefore, we can make use of the solution for the residual vectors from Equation (45) and insert them

in Equation (69) to obtain

σ2

0=vT

xk2+vT

yk1−avT

zk1−bvT

zk2

r. (70)

The variance-covariance matrix of the unknown parameters can be then derived by

Σˆ

xˆ

x=ˆ

σ2

0Qˆ

xˆ

x=





σ2

aσab σay0σax0

σba σ2

bσby0σbx0

σy0aσy0bσ2

y0σy0x0

σx0aσx0bσx0y0σ2







. (71)

To derive the variance-covariance matrix of the normalized vector components (64), we can

explicitly write the equations

an=ˆ

pˆ

a2+ˆ

b2+c2,

bn=ˆ

pˆ

a2+ˆ

b2+c2,

cn=c

pˆ

a2+ˆ

b2+c2.

(72)

Mathematics 2020,8, 1450 13 of 19

with

1. Following the standard procedure of the variance-covariance propagation in nonlinear

cases, we can write the Jacobian matrix







∂an

∂a

∂an

∂b

∂bn

∂a

∂bn

∂b

∂cn

∂a

∂cn

∂b







. (73)

Taking into account the variances and covariances of the line parameters ˆ

aand ˆ

bfrom (71)

Σˆ

aˆ

b=Σˆ

xˆ

x(1 : 2, 1 : 2) = "σ2

aσab

σba σ2

b#, (74)

we can compute the variance-covariance matrix of the normalized components

Σanbncn=FΣˆ

aˆ

bFT=





σ2

anσanbnσancn

σbnanσ2

bnσbncn

σcnanσcnbnσ2





. (75)

3. Numerical Examples

In this section we present the solutions for fitting a straight line to 3D point data using the TLS

approach from Section 2.1 and the newly developed WTLS approaches from Sections 2.2 and 2.3.

The utilized dataset for this investigaiton consists of

50 points and originates from the work of

Petras and Podlubny [

]. It has been utilized also by Snow and Schaffrin, see Table A1 in [

], for

a solution within the GH model, which will be used here for validating the results of the presented

WTLS solutions. Three different stochastic models will be imposed in the following:

equal weights, i.e., coordinate components

as equally weighted and

uncorrelated observations,

pointwise weights, i.e., coordinate components with same precision for each point and

without correlations,

general weights, i.e., correlated coordinate components of individual precision including singular

dispersion matrices.

3.1. Equal Weights

For the first case under investigation, we consider all coordinate components

as equally

weighted and uncorrelated observations, yielding the weights as shown in (3). The least squares

solution of this problem within the GH model, presented by Snow and Schaffrin [

], is listed in Table 1.

Mathematics 2020,8, 1450 14 of 19

Table 1. Least squares solution within the Gauss–Helmert (GH) model of Snow and Schaffrin [6].

Line Orientation Components Standard Deviation

Parameter ˆ

by=an0.219309 0.077523

Parameter ˆ

bx=bn0.677404 0.058450

Parameter ˆ

bz=cn0.702159 0.056575

Coordinates of a point on the line Standard deviation

Parameter ˆ

ay=y00.047785 0.121017

Parameter ˆ

ax=x0−0.067111 0.091456

Parameter ˆ

az=z00.049820 0.088503

A posteriori variance factor ˆ

σ2

00.7642885

A direct TLS solution for this problem can be derived using the approach presented in Section 2.1.

The results are shown in Table 2. Numerically equal results have been derived by using the direct

WTLS approach of Section 2.2 by setting all weights equal to one.

Table 2. Direct TLS solution from Section 2.1.

Line Orientation Components Standard Deviation

Parameter an0.219308632730 0.07752314583

Parameter bn0.677404488809 0.05844978733

Parameter cn0.702158730025 0.05657536189

A posteriori variance factor ˆ

σ2

00.76428828602

Comparing the solution with the one presented in Table 1, it can be concluded that the numerical

results for the parameters coincide within the specified decimal places. Regarding the small difference

in the numerical value for the variance factor, it is to be noted that the value in Table 2was confirmed

by two independent computations.

Furthermore, a point on the line can be easily computed using the equations of the functional

model (2), as long as the direction vector parallel to the requested line is known. Alternatively,

all the adjusted

points will lie on the requested straight line, which can be simply computed by adding

the computed residuals to the measured coordinates.

3.2. Pointwise Weights

For the second weighted case under investigation, we consider the coordinate components

of each point

to be uncorrelated and of equal precision. From the standard deviations listed in

Table 3, the corresponding pointwise weights can be obtained from (11).

Mathematics 2020,8, 1450 15 of 19

Table 3. Pointwise precision σxi=σyi=σzi=σifor each point Pi.

Point iσiPoint iσi

1 0.802 26 0.792

2 0.795 27 0.799

3 0.807 28 0.801

4 0.770 29 0.807

5 0.808 30 0.798

6 0.799 31 0.796

7 0.794 32 0.792

8 0.808 33 0.806

9 0.807 34 0.805

10 0.800 35 0.801

11 0.789 36 0.808

12 0.798 37 0.778

13 0.808 38 0.795

14 0.803 39 0.794

15 0.804 40 0.803

16 0.808 41 0.772

17 0.806 42 0.791

18 0.806 43 0.806

19 0.807 44 0.804

20 0.806 45 0.807

21 0.804 46 0.803

22 0.808 47 0.808

23 0.805 48 0.801

24 0.801 49 0.805

25 0.801 50 0.779

A direct WTLS solution is derived, following the approach presented in Section 2.2.

The determinant (23)

(d1−k)d4d5

d4(d2−k)d6

d5d6(d3−k)

=0,

can be built with the components

d1=213.505250528675,

d2=206.458905097029,

d3=198.273122927545,

and

d4=−20.9837621443375,

d5=−12.1835697465792,

d6=−81.6787394185243.

(76)

which leads to the solutions for the Lagrange multiplier

k=









115.0596477492 (ˆ

kmin)

218.3490470615

284.8285837425

(77)

The direct WTLS solution for the line orientation components is shown in Table 4.

Mathematics 2020,8, 1450 16 of 19

Table 4. Direct WTLS solution from Section 2.2.

Line Orientation Components Standard Deviation

Parameter an0.230818543507 0.07646636344

Parameter bn0.677278360907 0.05781967335

Parameter cn0.698582007942 0.05623243170

A posteriori variance factor ˆ

σ2

01.19853799739

The presented results are numerically equal to the iterative WTLS solution using the algorithmic

approach of Section 2.3, as well as the solution within the GH model.

3.3. General Weights

For the last weighted case in this investigation, we impose the most general case, i.e., correlated

coordinate components with individual precision resulting in a singular dispersion matrix. To obtain

such a matrix for our numerical investigations, we firstly solved the adjustment problem within the GH

model with an identity matrix as dispersion matrix. From the resulting 150

150 dispersion matrix of

the residuals

Qvv =





QvxvxQvxvyQvxvz

QvyvxQvyvyQvyvz

QvzvxQvzvyQvzvz





, (78)

computed as e.g., presented by Malissiovas ([

], p. 46), we take the variances and covariances

between the individual point coordinates, but not among the points, i.e., the diagonal elements of each

sub-matrix in (78) and arrange them in a new 150 ×150 matrix

QLL =





diag (Qvxvx)diag Qvxvydiag (Qvxvz)

diag Qvyvxdiag Qvyvydiag Qvyvz

diag (Qvzvx)diag Qvzvydiag (Qvzvz)







, (79)

with “

diag()

” denoting that only the diagonal elements are taken into account. This is an example of

pointwise variances and covariances, described as case (i) in Section 2.3, but now yielding a singular

dispersion matrix for the observations with

rank (QLL)=100.

Before deriving an iterative WTLS solution for this weighted case, we must check if

the criterion (66) for a unique solution of the adjustment problem is fulfilled. Therefore, we computed

the 100 ×100 matrix Wwith

rank (W)=100,

and the 100 ×4 matrix Awith

rank (A)=4.

The criterion ensures that a unique solution exists when using the presented singular dispersion

matrix, while

rank ([W|A]) =100, (80)

since

50 observed points are used in this example, cf. (66). As for all iterative procedures,

appropriate starting values for the unknowns must be provided. However, they can be obtained easily

by first generating a direct solution with a simplified stochastic model. The iterative WTLS solution for

the direction vector of the requested straight line is presented in Table 5.

Mathematics 2020,8, 1450 17 of 19

Table 5. Iterative WTLS solution from Section 2.3.

Line Orientation Components Standard Deviation

Parameter an0.225471114499 0.076563026291

Parameter bn0.677670055415 0.057791104518

Parameter cn0.699947192665 0.056127005073

A posteriori variance factor ˆ

σ2

00.798915322513

The presented WTLS solution has been found to be numerically equal to the least squares

solution within the GH model. Detailed numerical investigations of the convergence behaviour,

e.g., in comparison

to an adjustment within the GH model, are beyond the scope of this article.

However, in many numerical examples it could be observed that the iterative WTLS approach showed

a faster convergence rate compared to an adjustment within an iteratively linearized GH model.

4. Conclusions

For the problem of straight line fitting to 3D point data, two novel WTLS algorithms for two

individual weighting schemes have been presented in this study:

Direct WTLS solution for the case of pointwise weights, i.e., coordinate components with same

precision for each point and without correlations,

Iterative WTLS solution for the case of general weights, i.e., correlated coordinate components

of individual precision including singular dispersion matrices. This algorithm works without

linearizing the problem by Taylor series at any step of the solution process.

Both approaches are based on the work of Malissiovas [

], where similar algorithms have been

presented for adjustment problems that belong to the same class, i.e., nonlinear adjustments that can be

expressed within the EIV model. The approach presented in Section 2.1 provides a direct TLS solution

assuming equally weighted and uncorrelated coordinate components. The fact that this assumption

is inappropriate, e.g. for the evaluation of laser scanning data, has often been accepted in the past to

provide a direct solution for large data sets. With the newly developed approach in Section 2.2 it is

now possible to compute a direct WTLS solution at least for a more realistic stochastic model, namely

pointwise weighting schemes.

If more general weight matrices must be taken into account in the stochastic model, including

correlations or singular dispersion matrices, the presented algorithm of Section 2.3 can be utilized

for an iterative solution without linearizing the problem by Taylor series at any step, following

the algorithmic idea of WTLS. A criterion that ensures a unique solution of the problem when

employing singular dispersion matrices has also been presented, which is based on the original

ideas of Neitzel and Schaffrin [26,27], for a solution within the GH model.

Numerical examples have been presented in Section 3for testing the presented WTLS algorithms.

The utilized dataset of the observed 3D point data has also been employed by

Snow and Schaffrin [6]

for a solution of the problem within the GH model and originates from the study of Petras

and Podlubny [

]. The results of the presented algorithms have been compared in all cases with

existing solutions or the solutions coming from existing algorithms and have been found to be

numerically equal.

Author Contributions:

Conceptualization, G.M., F.N. and S.P.; methodology, G.M. and S.P.; software, G.M. and

S.W.; validation, F.N. and S.P.; formal analysis, G.M.; investigation, G.M. and S.W.; data curation, G.M. and

S.W.; writing—original draft preparation, G.M. and F.N.; writing—review and editing, G.M., F.N., S.P. and S.W.;

supervision, F.N. All authors have read and agreed to the published version of the manuscript.

Funding: This research received no external funding.

Acknowledgments:

The authors acknowledge support by the German Research Foundation and the Open Access

Publication Fund of TU Berlin.

Conflicts of Interest: The authors declare no conflict of interest.

Mathematics 2020,8, 1450 18 of 19

References

1. Joviˇci´c, D.; Lapaine, M.; Petrovi´c, S. Prilago ¯

davanje pravca skupu toˇcaka prostora [Fitting a straight line to

a set of points in space]. Geod. List. 1982,36, 260–266. (In Croatian)

Kahn, P.C. Simple methods for computing the least squares line in three dimensions. Comput. Chem.

1989,13, 191–195. [CrossRef]

Drixler, E. Analyse der Form und Lage von Objekten im Raum [Analysis of the Form and Position of Objects in Space];

Deutsche Geodätische Kommission bei der Bayerischen Akademie der Wissenschaften (DGK): München,

Germany, 1993; Volume C, No. 409. (In German)

Späth, H. Zur numerischen Berechnung der Trägheitsgeraden und der Trägheitsebene [Numerical calculation

of the straight line of inertia and the plane of inertia]. avn Allg. Vermess. Nachr.

2004

,111, 273–275.

(In German)

Kupferer, S. Verschiedene Ansätze zur Schätzung einer ausgleichenden Raumgeraden [Different approaches

for estimating an adjusted spatial line]. avn Allg. Vermess.-Nachr. 2004,111, 162–170. (In German)

Snow, K.; Schaffrin, B. Line Fitting in Euclidian 3D-Space. Stud. Geophys. Geod.

2016

,60, 210–227. [CrossRef]

Ghilani, C. Adjustment Computations, Spatial Data Analysis, 6th ed.; John Wiley and Sons, Inc.:

Hoboken, NJ, USA, 2018.

Niemeier, W. Ausgleichungsrechnung [Adjustment Computations], 2nd ed.; Walter de Gruyter: New York, NY,

USA, 2008. (In German)

Perovi´c, G. Least Squares (Monograph); Faculty of Civil Engineering, University of Belgrade: Belgrade,

Serbia, 2005.

10.

Pope, A. Some pitfalls to be avoided in the iterative adjustment of nonlinear problems. In Proceedings of the

38th Annual Meeting of the American Society of Photogrammetry, Washington, DC, USA, 12–17 March 1972;

pp. 449–477.

11.

Lenzmann, L.; Lenzmann, E. Strenge Auswertung des nichtlinearen Gauss-Helmert-Modells [Rigorous

solution of the nonlinear Gauss-Helmert-Model]. avn Allg. Vermess. Nachr. 2004,111, 68–73. (In German)

12.

Golub, G.; Van Loan, C. An analysis of the total least squares problem. SIAM

1980

,17, 883–893. [CrossRef]

13.

Van Huffel, S.; Vandewalle, J. Algebraic Connections Between the Least Squares and Total Least Squares

Problems. Numer. Math. 1989,55, 431–449. [CrossRef]

14.

Van Huffel, S.; Vandewalle, J. The Total Least Squares Problem: Computational Aspects and Analysis; SIAM:

Philadelphia, PA, USA, 1991.

15.

Malissiovas, G.; Neitzel, F.; Petrovic, S. Götterdämmerung over total least squares. J. Geod. Sci.

2016

,6, 43–60.

[CrossRef]

16.

Schaffrin, B.; Wieser, A. On weighted total least-squares adjustment for linear regression. J. Geod.

2008,82, 415–421. [CrossRef]

17.

Shen, Y.; Li, B.; Chen, Y. An iterative solution of weighted total least-squares adjustment. J. Geod.

2011,85, 229–238. [CrossRef]

18.

Amiri-Simkooei, A.; Jazaeri, S. Weighted total least squares formulated by standard least squares theory.

J. Geod. Sci. 2012,2, 113–124. [CrossRef]

19.

Snow, K. Topics in Total Least-Squares within the Errors-In-Variables Model: Singular Cofactor Matrices

and Prior Information. Ph.D. Thesis, Division of Geodetic Science, School of Earth Sciences, The Ohio State

University, Columbus, OH, USA, 2012; Report. No. 502.

20.

Malissiovas, G. New Nonlinear Adjustment Approaches for Applications in Geodesy and Related Fields; Ausschuss

Geodäsie der Bayerischen Akademie der Wissenschaften (DGK): München, Germany, 2019; Volume C,

No. 841.

21.

Guo, C.; Peng, J.; Li, C. Total least squares algorithms for fitting 3D straight lines. Int. J. Appl. Math.

Mach. Learn. 2017,6, 35–44. [CrossRef]

22.

Neitzel, F.; Petrovic, S. Total Least Squares (TLS) im Kontext der Ausgleichung nach kleinsten Quadraten

am Beispiel der ausgleichenden Geraden [Total Least Squares in the context of least squares adjustment of

a straight line]. zfv Z. für Geodäsie, Geoinf. und Landmanagement 2008,133, 141–148. (In German)

23.

Neitzel, F. Generalisation of Total Least-Squares on example of unweighted and weighted similarity

transformation. J. Geod. 2010,84, 751–762. [CrossRef]

Mathematics 2020,8, 1450 19 of 19

24.

Bronshtein, I.; Semendyayev, K.; Musiol, G.; Muehlig, H. Handbook of Mathematics, 5th ed.; Springer:

Berlin/Heidelberg, Germnay; New York, NY, USA, 2007.

25.

Boroviˇcka, J.; Spurný, P.; Keclíková, J. A new positional astrometric method for all-sky cameras.

Astron. Astrophys. Suppl. Ser. 1995,112, 173–178.

26.

Neitzel, F.; Schaffrin, B. On the Gauss-Helmert model with a singular dispersion matrix where BQ is of

smaller rank than B. J. Comput. Appl. Math. 2016,291, 458–467. [CrossRef]

27.

Neitzel, F.; Schaffrin, B. Adjusting a 2D Helmert transformation within a Gauss-Helmert model with

a singular dispersion matrix where BQ is of smaller rank than B. Acta Geod. Geophys. Montan. Hung.

2017,52, 479–496. [CrossRef]

28.

Wells, D.; Krakiwsky, E.J. The Method of Least Squares; Lecture Notes 18; Department of Geodesy and

Geomatics Engineering, University of New Brunswick: Fredericton, NB, Canada, 1971.

29.

Mikhail, E.M. Observations and Least Squares; Harper & Row Publishers: New York, NY, USA; Hagerstown,

MD, USA; San Francisco, CA, USA; London, UK, 1976.

30.

Petras, I.; Podlubny, I. State space description of national economies: The V4 countries. Comput. Stat.

Data Anal. 2007,52, 1223–1233. [CrossRef]

2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access

article distributed under the terms and conditions of the Creative Commons Attribution

(CC BY) license (http://creativecommons.org/licenses/by/4.0/).