ON THE ESTIMATION OF RIGID-BODY TRANSFORMATION FOR TLS REGISTRATION

The paper presents some considerations on the estim ation of rigid-body transformation which is the bas is for registration and georeferencing in terrestrial laser scanning. Two main issues are afforded, which in the opinion of the a uthor have been poorly investigated up until now. The first one concerns t he analysis of the stochastic model of Least Square s estimation of the transformation. In standard practice, the control p oints used for parameter estimation are considered as equally precise. Here two different approaches for weighting observations in laser scan and ground control points are proposed a nd evaluated through a set of Monte Carlo simulations. The second aspect refers to the reliability analysis. In particular the estima tion of a joint-adjustment including both geo-referencing parameters and geode tic n twork observations for the determination of g round control points is discussed and evaluated.


INTRODUCTION 1.1 Terrestrial Laser Scanning today
Terrestrial laser scanning (TLS) is coming out from its pioneer epoch and nowadays is applied to many practical fields.In many of them it plays a paramount role today.On the other hand, in the first decade since its appearance, many tests and experimentations have been carried out in order to assess instruments and methods, to refine processes, and to explore new applications for engineering surveying, deformation measurements, topographic works, geology, architecture, archaeology and cultural heritage documentation, 3D visualization, and Virtual Reality.On the other hand, some open problems still exist.These could be summarized in the following items: • standardization of instrument metrological features and methods for their evaluation; • definition of standard procedures for geo-referencing, scan planning, automatic co-registration of multiple scans; and • definition of application domains where the use of TLS is really worth, and where other techniques can be more fruitful (e.g.photogrammetry).
In the meantime, manufacturers have pushed up the diffusion of laser scanners, and purchasers have become familiar to the use of this technique for production work and research activity.
Very often the way they used TLS are limited to procedures proposed in vendor reference guidelines, without any critical analysis.This fact should be considered with great care, considering the recent increase of the interest for high precision applications, such as modelling of industrial sites (see Rabbani, 2006) and deformation measurement (see Vosselman and Maas, 2010).Here any step of the acquisition and processing of TLS data must be operated with the highest carefulness.
The basic geometric model implemented in laser scanning processing is the rigid-body transformation.A lot of work has been carried out so far on the estimation of the 6 parameters that are incorporated into this model (3 shifts and 3 rotations in space).Several aspects have been investigated, such as: evaluation of approximate values through linear or closed-form solutions, robust estimation, automatic labeling of corresponding points, integration of laser scans and images for mutual registration purpose.Such rigid-body transformation has been integrated by additional parameters to compensate for calibration errors (see Lichti, 2010).These solutions can be exploited to refine laser scanning measurements in applications where high accuracy is needed.Recently, different papers referred to a deeper analysis on the contribution of different error sources to the whole error budget of TLS measurement (Gordon and Lichti, 2005;Scaioni, 2005).Indeed, inclination of incident laser beam, diverse material reflections and other effects can result in discrepancies from the standard geometric model.The paper mainly affords two main issues, which in the opinion of the author have been poorly investigated up until now.
The first one concerns the Least Squares estimation of the rigidbody transformation.In the standard practice, points used for parameter estimation are considered as equally precise.On the other hand, if a scan is registered to ground or to another scan taken as reference, only points in one scan are usually considered as stochastic observations.To consider both problems, a formulation of the Least Squares rigid-body transformation estimation including a pointwise weighting in both reference systems (scan to ground or scan 1 to scan 2) is applied in Section 2. Some considerations based on different point configurations are reported to show when this approach is worth to be used.
The second aspect refers to the reliability analysis, i.e. the chance to detect gross errors in observations.Which are the factors that can improve this important property?In many closerange applications, ground control points (GCP) are measured within a local geodetic network.Errors in geodetic and laser scanning measurements can be at the same order.Consequently, a joint-adjustment of both sets of observations would be an option to increase the inner reliability.Some theoretical considerations about this solution and a simulated case are proposed in Subsection 3.4.

Definitions
In this paper the problem of two or more scan registration is focused.The term registration means the computation of the transformation(s) allowing all scans to be referred to the same reference system (RS).In the following, the instrumental RS of the scan i-th is termed Intrinsic RS (IRS i ).If an external geodetic datum is established, this is referred to as Ground RS (GRS); in this case, the registration of one or more scans into the GRS is addressed to as the geo-referencing problem.
The geometric model usually adopted for registration or georeferencing is a 3D rigid-body transformation, possibly integrated by a set of parameters for laser scanner calibration: where: x i : vector of 3D coordinates of a point in IRS i ; x j : vector of 3D coordinates of a point in IRS j (or in GRS); R ij : 3×3 rotation matrix parameterized through 3 rotation angles; t ij : 3D shift vector; ∆x i AP : vector including additional parameters for calibration in IRS i ; ∆x j AP : vector including additional parameter for calibration of points in IRS j (∆x j AP =0 if IRS j≡ GRS).
The computation of unknown parameters is carried out on the basis of common features (usually points are used) in both RSs.
On the other hand, in practical applications, the geometric model ( 1) is rarely integrated with further additional parameters to correct errors due to non-modelled systematic effects (∆x AP ).However, when this happens (see e.g.Lichti, 2008), the calibration is performed at a preliminary stage by measuring a test-field made up of tens of GCPs or other kinds of features (e.g.planes), and the results applied to correct the x 2 i vector in (1).The estimation of the unknown parameters of geometric model (1) is usually carried out by Ordinary Least Squares (OLS), after linearization (Wolf and Ghilani, 1997).This requires the computation of approximations for all parameters, which can be derived by a twofold strategy.The first one relies on the use of an algorithm which can perform the solution of Eq. (1) without any initial values.Three main approaches are adopted based on Procrustes analysis (Beinat et al., 2000), Hamilton's Quaternions (Horn, 1987), Orthonormal transformations (Horn et al., 1998).At the first stage, often the automatic search for correspondent features in both IRS i and IRS j is applied (Bornaz et al., 2003).Furthermore, a high breakdown-point estimator (RANSAC or Least Median Squares, for instance - Rousseuw and Leroy, 1987) is often implemented to cope with gross errors.Finally the estimate with OLS is applied by using the linearized model ( 1).An automatic testing or manual check of remaining small outliers can be included at this stage.On the other hand, when a TLS is used along with direct geo-referencing technique (Scaioni, 2005), approximations of parameters are already known.

Functional model for Ordinary Least Squares
The functional model for OLS defines the linear (or linearized) relation between observations and parameters to estimate.According to the formulation based on observation equations, points measured in scan i-th (x i ) are organized in the vector y 0 .These are function of the points measured in scan j-th (x j ) as well as the corrections to registration parameters (dβ 1 ), to be estimated through eq. ( 1) after linearization.To consider also points in scan j-th as stochastic variables, they are introduced as further unknowns (dβ 2 ).Consequently two groups of equations are setup.The first comprehends the observation equations: being A 1 and A 2 the design matrices corresponding to registration parameters (dβ 1 ) and to points in scan 2 (dβ 2 ); v is the vector of residuals which are minimized in OLS; c is the known term vector.To introduce uncertainty also on points in scan j-th, a second group of pseudo-observation equations should be included in the functional model: To account for the uncertainty of observations, two weight matrices are adopted (W i , W j ) for points observed in scans i-th and j-th, respectively.Three relevant options for the weight matrices are thoroughly analysed in Paragraph 2.1.2.

Stochastic model
In the most applications, only points in the scan i-th are considered as stochastic variables (cases 'a' and 'b').
In the case 'a', all measured points have the same precision.The standard model for error distribution is the hysotropic Gaussian distribution, assumption that only approximately corresponds to reality.Indeed, the use of the same weights for all points is not correct.Indeed, measurement errors in laser observations (angles and ranges) are not hysotropic.Moreover, some error sources in laser measurements are not normally distributed (e.g. the effect of laser beamwidth and surface reflectivity - Lichti et al., 2005).Although the Central Theorem of Statistics supports the use of the Gaussian model, this is only an approximation.Indeed, the most error sources are neglected in the stochastic model because they are difficult to quantify.Some results have been published about the effects of surface orientation with respect to laser beam (Bae et al., 2009).However, this approach is difficult to follow in the practice.Here a simple solution for weighting points in a scan is used (case 'b'), by adopting the 'positional' model proposed in Alba et al. (2010).The precision of each point k in the scan can be expressed through the covariance matrix: The matrix C k has been constructed by using the standard deviations of the laser scanner measurements (σ ρ ,σ α ,σ θ ), including range ρ and both horizontal (α) and vertical angles International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XXXIX-B5, 2012 XXII ISPRS Congress, 25 August -01 September 2012, Melbourne, Australia (θ).When an algorithm for automatic target measurement is adopted, a further contribute to the error budget can be added up by introducing related variances and covariances (indicated with sub-index 'FM').The matrix (4) in polar coordinates has to be propagated to obtain the one in Cartesian coordinates (C′ k ) needed for the estimation of parameters of Eq. ( 1).Then the covariance matrix (C) of all the observed points in the scan will result from the composition of all matrices C′ k into a blockdiagonal structure.Consequently the weight matrix will be W= C -1 .On the other hand, other more involved error models can be used, like that proposed in Bae et al. (2009).In any case, when operating with real datasets, the effectiveness of the adopted method depends upon how well the error model fits real observations.On the other hand, simulated data allows to understand which contribute can be achieved from the use of weighted observations.A third case ('c') takes into consideration also weighted observations in both scans i-th and j-th.Both weight matrices W 1 and W 2 can be defined as in case 'b'.This formulation is adequate to describe the co-registration of two laser scans in the RS of one of them.If the RS of one scan is given by GCPs, these can be weighted by adopting the estimated precision obtained from the LS adjustment of the geodetic network.

Monte Carlo analysis of stochastic models
Five simulated datasets have been setup with various numbers of points and spatial distribution (Fig. 1).In each of them, a laser scan has been acquired from a central position with respect to a set of targets to be used as GCP for geo-referencing.
Instrument precision is that typical of a medium range phaseshift modern scanner (σ ρ =±5 mm; σ α =σ θ =±0.05 mrad).The 'positional' error model ( 4) has been adopted to simulate measurement errors.The precision of target measurements or the presence of systematic errors has been neglected here.In the cases 'a' and 'b' the GCPs have been considered as error-free.
In the case 'c' they have been weighted by using the estimated theoretical accuracies from geodetic network adjustment.
A Monte Carlo simulation (Robert and Casella, 2004) has been applied to all datasets to compare diverse weighting strategies.The evaluation of results is possible because the true value of points is known here.This makes possible to analyse either residuals on error-free independent check points and parameters.This is obviously a great advantage to verify the adequacy of different stochastic models under investigation.On the other hand, simulations do not allow one to assess how the assumed models fit with real data.Monte Carlo simulation is based on the repetition of a high number of trials where input data are randomly extracted based on a probabilistic distribution.Here the simulated observations consist on a set of corresponding points which have been measured in both IRS and GRS.A set of preliminary experiments have been carried out to setup the minimum number of trials (n tr ) giving enough stability to the results.Once a value n tr =10000 is decided, 10 sets of experiments with the same number of trials has been repeated to assess the repeatability of results.This assumption has been verified, considering that all outcomes differ less than the assumed data uncertainty.

Results of simulations
Configurations from A1 to A4 have been designed to progressively reduce stability.This result has been obtained by either decreasing the total number of points, and by weakening their spatial distribution.In configuration A5 a small set of points (#12) like in A4 has been setup, but with a stronger geometry.
The quality of geo-referencing has been evaluated by computing residuals on error-free check points.Two groups of check points have been analysed.The first collects all points of dataset A1 and allows evaluating the global quality of the estimated geo-referencing.The second group comprehends only those points contained in the region internal to the GCPs.Indeed, best practices always suggest selecting GCPs around the area to survey.
Figure 1.-GCP configurations adopted in Monte Carlo simulations to assess different stochastic models.In general, computed residuals on check points have resulted smaller in the second subset than in the first group, as expected.
Although there is not full evidence that a stochastic model is absolutely better than another, two conclusions can be drawn.In the most stable configurations (A1-A2) all stochastic models worked well.As far as the configuration becomes weaker, the case 'c' showed slightly better results than 'b.'

Criterion matrix
The simulations' outcomes suggest checking the spatial distribution of points to detect possible weak configurations.Such kind of criteria is also important in estimation techniques (e.g.RANSAC or Least Median of Squares) based on random extraction of several minimal point datasets to be used for registration.Such assessment can be based on the comparison of the estimated covariance matrix of parameters (C xx ) with an ideal one called criterion matrix H (Baarda, 1973).This test allows one to establish if the current point configuration is better than expected, according to H.Here the approach proposed in Förstner (2001) has been modified to account for the different number of points in the real (n C ) and the reference case (n H ): The comparison is accepted depending upon the largest eigenvalue of matrix K: Both expressions in Eqs. ( 5) and ( 6) can be easily worked out when the number of parameters is small (6 in the case under study).The introduction of the weighting coefficient n C /n H is motivated by the observation that if H is setup based on an ideal configuration, the number of points in here might strongly affect the results with respect to their spatial distribution.Alternatively, Förstner (2001) suggests to compute H based on some theoretical considerations on variances and correlations between parameters.
In Table 2 the results of test (6) for the 5 GCP configurations in Figure 1 are shown.In addition, some subsets of points have been derived from the largest dataset A1.While a homogenous spatial distribution of points has been kept, their number has been progressively reduced.

case #points
-Results of comparison of different point configurations with a criterion matrix from case A1.

Basics concepts
Reliability refers to the chance to detect gross errors in the target coordinates measured in the laser scan and adopted for geo-referencing purpose.The same concept could be extended to include also errors in GCP coordinates, but this case is not considered here.
The problem is to evaluate how big can be the error on estimated geo-referencing parameters of a scan, according to the largest theoretical error in observations that cannot be theoretically detected by data snooping (Baarda, 1968).The inner reliability relates the expectation E(|v i |)≠0 (also called non-centrality parameter δ 0 ) of the normalised correction corresponding to a gross error ∆l i which is locatable with a given probability 1-β: International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XXXIX-B5, 2012 XXII ISPRS Congress, 25 August -01 September 2012, Melbourne, Australia where r i is the local redundancy of the observation l i , and σ li the corresponding mean square error.If the gross error is smaller than the inner reliability E(|∆l i |), the probability that it will not be detected in a test with normalised corrections is high.Consequently, this error will result in a deviation of the estimated parameters x from LS adjustment.The corresponding outer reliability is then given by the expression: ( ) ( ) Vector E(|∆l|) is zero except on the element corresponding to the observation affected by gross error ∆l i .

Analysis of simulated configurations
A reliability analysis of the configurations adopted in Section 2 has been carried out.A non-centrality parameter δ 0 = E(|v i |)=4.0, a significance α=0.01 and a power β=0.93 have been setup.The corresponding rejection threshold for data snooping is k=2.56.
In Table 3 some results on the inner reliabilities computed for points in configurations A1, A2, A3, A4, and A5 are reported.
As can be seen, in the less stable configurations (A4-A3-A5) the largest inner reliabilities can reach some centimetres.
Assuming that in the dataset is present only one gross error which is equal to the largest inner reliability value, this will lead to a bias in the estimated geo-referencing parameters.This can be computed based on Eq. ( 8).To better understand the effect of this error on the final point coordinates, each biased set of parameters has been used to compute the discrepancies on the check point coordinates.As can be seen in Table 3, the effect on check points of the largest undetectable error in observations is quite relevant for weak configurations.On the other hand, no significant effects can be noticed in the strongest ones (A1-A2).This outcome highlights the need of major care on reliability analysis during geo-referencing, especially in high-precision applications.

Leverage points
Leverage points in regression estimation are meant to be points which have a strong influence on the estimated result.This is due to their geometric position, especially in the case of nonhomogenous or unbalanced data.Such points are characterized by a low local redundancy, and the corresponding residuals can mask larger errors.The corresponding inner reliability will be larger, according to Eq. ( 7).Table 3 reports the inner reliabilities computed for all datasets in Section 2. It is evident how no leverage points appear in configurations A1 and A2, due to the large number of data and their good spatial distribution which result in quite even inner reliabilities.The result is different in smaller and asymmetric configurations like A3, A4 and A5.Here there are some points with higher inner reliabilities than average values.Two new small datasets (see Fig. 2) have been setup for the purpose of investigating the existence of leverage points (LEV1 and LEV2).The lower rows of Table 3 show how dramatic is the problem here.If in general the average inner reliabilities are large due to the poor data redundancy, some observations lead to undetectable errors up to several centimetres, as can be seen in the rightmost three columns of Table 3.Here the maximum residuals on check points depending on the error on the observation with the largest inner reliability have been computed based on Eq. ( 8).In reality, the situation can be also worse because observations in TLS are strongly correlated among them.Indeed, here the coordinates of targets in a scan have been considered as independently observed data, but in fact they are not.

Joint-adjustment of laser scan and geodetic data
A further chance to increase the local redundancy of laser scan observations is to combine the adjustment for computing georeferencing parameters to the one of the geodetic network for the determination of GCP coordinates.The procedure is similar to joint-adjustment that was experienced in aerial photogrammetry.
The design matrix A is made up of different sub-matrices.Observation equations coming from linearization of Eq. ( 1) contribute to both sub-matrix A p (size 3n×6, where n is the number of points), containing the coefficients of georeferencing parameters, and to A GCP (size 3n×3n), with the coefficients of GCPs.Three different kinds of theodolite measurements can be adopted in a geodetic network for working out GCP coordinates: azimuth and zenith angles, slope distances.The geodetic datum can be setup thanks to the knowledge of the stations' coordinates, or by arbitrarily fix them to establish a minimum constraint.This second group of observation equations gives rise to sub-matrix A g (size h×3n, where h is the number of geodetic observations).The full design matrix A and the weight matrix W of joint-adjustment is then given by: where sub-matrices W TLS and W g correspond to laser scanning and geodetic observations, respectfully.The redundancy matrix can be computed as follows: Table 4. -Features of independent and joint-adjustment for estimating geo-referencing parameters and GCP coordinates.

CONCLUSIONS AND FUTURE DEVELOPMENTS
In recent years terrestrial laser scanning (TLS) techniques have achieved a great success in many domains.For the most this is due to the straightforward acquisition of 3D data and to the simplification of the registration procedures with respect to photogrammetry.In many practical applications, today TLS allows users to accomplish projects that in the past required a much higher effort.On the other hand, a different approach should be followed in the case of high-precision applications, likewise deformation monitoring.
Although much work has been carried out on modelling systematic errors, automation of laser scan registration, and point cloud segmentation, less attention has been put on design of data acquisition and quality assessment.In close-range photogrammetry, paramount work on these topics was carried out in the '80s and '90s.This paper is a short attempt to draw the attention of TLS researchers and practitioners on some of these concepts.The development of stochastic models that can be suitable to deal with real data is an important issue, especially when dealing with weak configurations in geo-referencing.The reliability analysis has been demonstrated to play a fundamental role in high-precision applications.Here the case of single scan registration to ground has been considered.The same analysis could be extended to the case of several scans to be registered in a common frame.

Figure 2 .
Figure 2. -GCP configurations with leverage points; red triangles are GCP, blue circle check points that are outside the bounding box of GCP to simulate a typical scenario in deformation measurements.
Table 1.-Results of Monte Carlo simulations of a scan geo-referencing based on different GCP configurations and stochastic models ('a': uniform weighting of laser scanner observations only; 'b': weighting laser scanner observations based on the 'positional' model; 'c': using also weights for GCPs).Results refer to the use of two different groups of check points.

Table 3 .
-Results of reliability analysis.