THE IMPACT OF THE DISTRIBUTION OF OBSERVATIONS ON TERRESTRIAL LASER SCANNER SELF-CALIBRATION QUALITY

Terrestrial laser scanning (TLS) instruments find routine use for a range of precision engineering measurement applications. Depending on the accuracy requirements for a specific project, the instrument may require self-calibration to determine systematic error model terms. One of the goals of first-order network design for self-calibration is to acquire observations throughout the full instrumental field-of-view. Experience calibrating TLS instruments has demonstrated that while this goal can be achieved for horizontal deflection angle observations, it is seldom realized for the vertical angle observations. This paper presents results from a preliminary investigation into the impact of the distribution of vertical angle observations on the estimation of two critical systematic error parameters in TLS instruments: the collimation axis error and the trunnion axis error. First, a model to characterize the empirical observation distributions is developed. The model is a function of a single shape parameter that quantifies observation dispersion. Then, a means to estimate the impact of the distribution on the parameter estimation is developed. Results from six real datasets show the distribution model characterizes the overall general trend of the observations. Simulated results show the relative independence of the collimation axis error and the strong dependence of the trunnion axis error on the shape parameter.


Terrestrial Laser Scanner Calibration
Sensor modelling and calibration are recognized as important quality assurance processes for terrestrial laser scanning (TLS) instruments used for precise engineering measurement projects. Systematic errors inherent to TLS instruments must be identified and the model parameters estimated so that acquired point clouds can be corrected in order to maximize accuracy. Self-calibration has been identified as an effective procedure to determine these parameters. TLS self-calibration methods first emerged more than a decade and a half ago (Gielsdorf et al., 2004) and recent work has shown that they are still necessary (Lichti et al., 2019). Self-calibration can be performed in a laboratory setting (Lichti, 2007) or on site (Abbas et al., 2014) using either signalized targets (Muralikrishnan et al., 2015;Reshetyuk, 2010) or geometric primitives inherent to the scene (Chow et al., 2013;Medić et al., 2017).
A geometrically strong network is generally regarded as essential to the accurate estimation of systematic error parameters. Lichti (2007) reports that a large elevation angle range is needed for the estimation of certain error terms. Moreover, Lichti (2010) used simulation to investigate the impact of the range of elevation angle measurements on the precision of and the correlation between certain model variables. However, the exact impact of the distribution of these observations on the self-calibration solution is not known. This work focuses on modelling and quantifying the impact of the distribution of angular observations within a TLS self-calibration network.

Camera Calibration Analogy
To better illustrate this problem, an analogy to camera calibration is drawn. This is a well-documented process, e.g. (Luhmann et al., 2018), to determine the interior orientation parameters of a camera. Strong first-order network design measures are required to reduce the effects of projective compensation (Fraser, 1997) and to maximize parameter quality. One important design feature is to fill the image format or, equivalently, the camera's field-ofview (FoV), with image point observations. This aim of this measure is to improve the quality of the estimated radial lens distortion coefficients. The gradient of the distortion curve is typically greatest at the edges of the format, so observations near the periphery of the image plane have a strong impact on distortion parameter precision.
A simple example is illustrated in Figure 1. For the sake of clarity, the camera calibration is distilled to a curve fitting problem to determine a single parameter of the radial lens distortion profile, k1, from 2D point coordinates. Two sets of 100 observation points are shown. The location of each point was randomly drawn from a 2D uniform density function. Whereas the first set of observations spans the full image format (8 mm x 6 mm), the second set only covers 80% of the image plane.
The differences in error envelopes of the two curves clearly show the strong influence of the point distribution on parameter quality. A 20% reduction in data coverage leads to a 100% increase in the standard deviation of the distortion profile, δr.
(The exact figure depends on the randomly-simulated point distribution, but from many repeated trials the factor was found to be about 100%.)

Terrestrial Laser Scanner Observations
Returning to the laser scanning problem, TLS instruments collect range observations, ρ, by time-of-flight laser range finding in a spherical imaging geometry ( Figure 2). The emitted laser beam is deflected in nominally equal increments of elevation angle, α, within a vertical profile. The rotating instrument head allows the capture of a series of vertical profiles in nominally equallyspaced increments of horizontal direction, θ. The spherical observations for point i appearing in scan j can be parameterized in terms of the scanner space Cartesian coordinates (x, y, z) as follows: Note that the observations have been augmented with systematic error correction terms (∆ρ, ∆θ, and ∆α) that are described in the next sub-section. Different instrument architectures exist (Staiger, 2003). The focus here is on the panoramic geometry in which the respective angular fields-of-view for the horizontal direction, θ, and the elevation angle, α, are (0° ≤ θ < 180°) and (-αmin ≤ α ≤ 180°+αmin). The parameter αmin ranges from 60° to 70°, depending on the instrument. As a result of this architecture, observations are collected "in front of" and "behind" the instrument. Thus, observations at the horizon are collected at two locations: α=0° and α=180°.
As mentioned, these APs can be estimated by the target-based self-calibration method. Scans are captured from different locations within a controlled environment comprising an array of signalized targets. The location and orientation of the scans and the distribution of the targets governs the distribution of the observations. It is hypothesized that, as demonstrated in the camera calibration example, the distribution of observations within the instrument FoV governs the AP solution quality.

Terrestrial Laser Scanner Observation Distribution
Histograms of the angular observations from six TLS selfcalibration datasets captured with five different panoramic instruments are shown in Figures 3 and 4. They were acquired in different indoor 3D calibration target fields. These datasets are described in greater detail in Section 3. Each dataset comprises observations captured from multiple instrument locations within the respective target field.
Clearly, the horizontal direction observations (θ) are generally uniformly distributed throughout the horizontal FoV. However, despite the placement of many targets on the ceiling and the floor The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII- B1-2020, 2020XXIV ISPRS Congress (2020 in the calibration rooms, the elevation angle observations are not uniformly distributed. Instead, they are clustered in two nearlysymmetric lobes centred about the horizon (at α=0° and α=180°).
The data are very sparse at large elevation angles (e.g. α = -αmin, 90°, and 180+αmin). Two of the basic model APs, the collimation axis error, b1, and the trunnion axis error, b2, strongly depend on the elevation angle. They respectively vary with the secant and tangent of elevation angle. Given the data characteristics in Figure 4 and the steep gradient of the functions involved, one might expect a strong dependence of the precision of b1 and b2 on the distribution of the elevation angle observations.
With this background, the aims of this work are to: 1. Develop a model for the distribution of elevation angle observations from real datasets; 2. Examine how the distribution of the elevation angle observations influences the estimation these two APs; 3. Develop a means to predict the quality of two APs, the collimation axis error and the trunnion axis error, given the hypothesized model for the elevation angle distribution; 4. Validate the model prediction of AP quality by comparison with estimates from real dataset selfcalibrations; and 5. Propose an observation distribution and, by extension, a calibration target field design, that maximizes the precision of the two APs.
This paper is a report on preliminary efforts to address aims 1, 2 and 3.

Observation Distribution Modelling
As can be seen in Figure 4, a wide variety of empirical histogram shapes exists. Whereas some lobes are very compact, others exhibit greater dispersion. The compactness/dispersion indicates the range of elevation angle observations over which data were  The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B1-2020, 2020 XXIV ISPRS Congress (2020 edition) collected for the self-calibration. Gaps exist in the histograms due to target field design, sample size and the choice of bin width. Moreover, the distribution depends on the target field layout and TLS instrument location, as mentioned previously. This complicates the development of a model that characterizes the observation distribution.
A straightforward but realistic simulation has been performed to illustrate the role of target field layout and instrument location in determining elevation angle observation distribution. The real datasets were used to generate the 2D TLS calibration target field shown in Figure 5 so that the dimensions and the number of targets (200) are realistic. Two different sets of three TLS instrument (αmin=70°) locations are also shown. All targets were assumed to be visible from each instrument location. The targets and the instrument locations are confined to the xz plane, Spherical observations were derived from the Cartesian coordinates with y=0. The resulting histograms, each comprising 600 observations, are shown in Figure 6. Clearly, both exhibit similar general behaviour, but they differ considerably in terms of fine-scale details.
Density estimation is well-studied topic. In the univariate case, the goal is to determine the probability density function of a single random variable from a finite set of observations (Searle, 1987). In this work, the function should exhibit certain traits. First, it should capture the overall general behaviour of the elevation angle observation distribution. Second, it should be relatively simple in terms of its functional form and involve as few parameters as possible in order to keep the mathematics described in the next sub-section tractable. Ultimately, the choice represents a trade-off between model complexity and goodness of fit to the data. Methods such as kernel density estimation can be used to obtain a very good fit to each of the distributions shown in Figures 4 or 6. This is not the aim here. The aim is to capture the overall trend under the hypothesis that a function describing the general behaviour will give insight into impact of the observation distribution on AP estimation quality.
A simple parametric model that captures the general distribution behaviour is therefore sought. Although many possibilities exist, the biomodal raised cosine density function has been investigated here: It is a function of a single parameter, α0, which represents the cutoff elevation angle for the observations and governs the shape of the distribution. This parameter also defines the extents of the distribution. The two lobes are centred at the horizon, 0° and 180°. Raised cosine distributions for several values of α0 are shown in Figure 7.
Estimation of α0 from sample data is straightforward. The variance, σ 2 , is computed with the mean of either 0° or 180°, Figure 5. Simulated 2D target field and instrument locations.

Figure 6. Elevation angle histograms for the 2D network for instrument location set 1 (top) and set 2 (bottom).
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B1-2020, 2020 XXIV ISPRS Congress (2020 edition) whichever is closer to a particular sample. The shape parameter is computed from the variance as follows

AP Quality Prediction Model
To develop a model for AP quality, the correction model for the collimation axis and trunnion axis errors can be investigated separately from other model variables. This independent analysis is possible because these two APs are effectively decoupled from the other parameters in the panoramic scanner self-calibration solution when the geometric network design is strong. From the six datasets used herein, the maximum correlation coefficient between the two APs in question and the exterior orientation parameters was 0.37. The multivariate estimation problem thus reduces to a simpler regression task.
The observation equation observation i is obtained by combining Equations 2 and 5 and rearranging and by making the following substitution Note that the scanner subscript j has been dropped since the analyses can be performed without reference to object space.
For a set of n observations, the system of observation equations is given by tan where A, x, y and v are the design matrix, parameter vector, observation vector and residual vector, respectively. The errors in the θ′ observations are assumed to be drawn from the same population, so the weight matrix P can be defined as a scalar matrix: The elevation angles, αi, are assumed to be error free.
The least-squares normal equations are given by The normal equations matrix, N, is of particular interest since its inverse is the cofactor matrix of the APs, which quantifies their quality. In this problem N has the following analytical form Under the assumption of uniform sampling within a vertical profile, ∆α is the increment between two successive observations. Multiplying the normal equations by this scalar results in The normal equations are further multiplied by the function p(α), which represents the distribution of the vertical angle observations.
As the number of observations becomes large, that is in the limit as n→∞, the summations above can be rewritten using the definition of the definite integral All integrands must be continuous on (a,b), where a=-90° and b=270°. This condition is not strictly met since the secant and tangent functions are infinite at α=90°. So, the two branches of each integral must be evaluated separately. Note that the offdiagonal term is equal to zero if the limits are symmetric since the product of the secant and the tangent functions reduces to the odd-symmetric sine function.
Equation 18 is an expression for the normal equations matrix that explicitly models the distribution of the observations. The inverse of this matrix is, of course, the cofactor matrix of the parameters, Qx, which contains parameter quality information.
Thus, the evaluation of the integrals and subsequent matrix inversion allows investigation of AP precision as a function of elevation angle distribution. The model developed using the procedure described in Section 2.1 can be substituted into Equation 18 to investigate the impact of the shape parameter, which models the limits of data acquired in a self-calibration.
Furthermore, other data distributions based on different target field designs can be investigated. Ultimately, this method can be used to identify an optimal distribution. This approach avoids the use of simulation techniques that necessitate performing many trials with artificially-generated observations. Instead, this completely analytical approach provides a more elegant means to estimate quality as a function of data distribution. It is therefore a powerful first-order network design tool.

EXPERIMENT DESCRIPTION
The methods developed in Section 2.1 have been applied to the six TLS self-calibration datasets depicted by the histograms in Figures 3 and 4. Some of the pertinent details are provided in Table 1. All were captured under controlled temperature and lighting conditions in a purpose-built, indoor 3D array of highcontrast targets. As can be deduced from the tabulated metadata, each dataset comprises multiple scans of many targets and is highly redundant. Further details about these can be found in related publications. Datasets 1 and 2 are described in (Lichti, 2007). Datasets 3 and 4 are documented in (Lichti et al., 2019). Datasets 5 and 6 are respectively described in (Al-Manasir and Lichti, 2015) and (Chow et al., 2013). For each dataset the shape parameter α0 was estimated and the quality of the model fit was evaluated. In the second experiment, the cofactor matrix described in Section 2.2 was formed and the integrals were evaluated as a function of distribution shape. The observation variance, σ 2 , was assumed to be unity.

Distribution Modelling
Estimates of the shape parameter are presented in Table 2. There is a wide range of values for the shape parameter α0: from 50° to more than 80°. The estimated distribution functions are superimposed on the histogram from Figure 4 in Figure 8. The raised cosine does capture the general shape of the histograms. However, there are large deviations at the horizon where more observations were captured than the raised cosine model can predict. Deviations are also visible where gaps occur in the histograms. These outcomes were anticipated from the simulation reported in Section 2.1.
The goodness-of-fit is graphically depicted in Figure 9, which shows the cumulative distribution functions (CDFs) determined empirically from the elevation angle observations and the raised cosine CDFs computed from the estimated shape parameter. Note that the data above the zenith have been mapped to (-90°<α<90°) for the sake of clarity. The unimodal and bimodal raised cosine functions are related in amplitude by a factor of 2 so there is no loss of generality by performing this transformation. Visual inspection shows that the model fit is better for some cases than others, but all follow the general trend. The Kolmogorov-Smirnov (K-S) test is a non-parametric test performed to examine the quality of fit of an empirical distribution to an hypothesized distribution. The test statistic is the maximum difference in cumulative probability between the two distributions (Kanji, 1999) Table 2 gives the results of the K-S test for all six datasets. None of them pass the test at the 5% level of significance. Although the aim here is to model the general behaviour rather than the fine details, more investigation into model choice may be warranted. Figure 10 shows the behaviour AP quality versus the shape parameter, α0. The square root of the cofactor, q, is used as an indicator of precision for each AP. The collimation axis error precision is only mildly affected by the distribution of observations. It is effectively constant up to about 48°, where it deviates by 5% or less from the precision at 5°. After that, the precision gradually improves by just under 22% at 85°, the maximum computed elevation angle.

AP Precision Prediction
The trunnion axis error precision is strongly affected by the observation distribution. The precision is very poor at low elevation angles where the tangent function is only weakly observable. The trunnion axis error is indeterminate in the limit as α0→0°. It is more than 30 times worse than that of the collimation axis error at 5°. It improves considerably, though, as the range of elevation angle observations increases. At 85° the trunnion axis error precision is only 1.6 times worse.

CONCLUSIONS AND FUTURE WORK
The quality of systematic error parameters estimated by instrument self-calibration depends on the distribution of observations throughout the FoV. The distribution of observations is a function of network design. In TLS selfcalibration networks, elevation angle observations tend to be clustered around the horizon, which may have an impact on some APs due to their strong gradients near zenith and nadir. In this work, a model for the distribution of TLS observations from prior calibration datasets has been developed to investigate how additional parameter quality is influenced by the distribution. The raised cosine function, which depends only on one shape parameter, was investigated for this purpose. With focus on the collimation axis error and trunnion axis error APs, a means to predict their precision from the least-squares normal equations was developed. The application of the hypothesized distribution function and quality prediction method revealed that the collimation axis error is insensitive to the shape parameter. The trunnion axis error, though, was found to be strongly dependent on the observation distribution. The latter finding suggests that more observations at the extents of the field of view are needed to improve quality. Although this follows conventional general wisdom, using this modelling approach will allow future work on the development of an optimal distribution function. This will  The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B1-2020, 2020 XXIV ISPRS Congress (2020 edition) provide informed guidance for TLS calibration network design and improve trunnion axis error estimation. Additional datasets will be acquired for further testing. Moreover, the model estimates for AP precision will be compared to actual results from real TLS self-calibration datasets.
The observation distribution model was chosen for mathematical simplicity and the ability to model the general observation trend. The selected model did not fit any of the datasets at the 95% confidence level. Future work could be devoted to investigating other, more detailed models. In addition, the sensitivity of AP estimates on model shape can be investigated.