PROBABILISTIC FEATURE MATCHING APPLIED TO PRIMITIVE BASED REGISTRATION OF TLS DATA

Many industrial applications require dense point clouds of the installations. Acquisition of the rooms, filled with many objects, of an industrial scene leads to many Terrestrial Laser Scanner (TLS) stations. A precise registration of all the per-station point clouds is crucial for the required accuracy of 1-2 cm of final data. Targets and tachometry, current best practice for registration, slows down the survey and limits the number of campaigns. Indoor geolocation system are faster but do not reach the final required accuracy. Otherwise, 3D primitives can be automatically extracted from the dense point clouds and possibly used for registration. In a four step primitive-based registration, Acquisition Reconstruction Matching Solving, the matching is crucial. This article presents a probabilistic test for 3D lines matching using a priori distributions of approximated transformations. The stochastic model of approximated transformations and resulting uncertain lines is introduced. A test is performed on a real dataset of an industrial scene and the results are analysed. Improvements of the presented test and matching framework are also discussed.


INTRODUCTION
Maintenance and revamping of large industrial installations commonly require 3D As-Built models of their geometries.The 3D reconstruction of these complex scenes can be achieved from an accurate and dense point cloud.Current laser scanners can deliver a high quality point cloud for each station.The correct registration of this Terrestrial Laser Scanning (TLS) data is crucial for the accuracy of the global point cloud.Accurate tachometric measurements based on widely distributed targets can provide the required accuracy of 1-2 cm.This protocol is now the bottleneck of wide survey campaign (over 1000 stations).
An alternative approach exploits geometric primitives as "natural targets", which are already present in the scene (points, lines and planes respectively stemmed from spheres, cylinders and walls).These primitives are automatically, or semi-automatically, reconstructed from per-station point clouds.The framework presented in (Hullo et al., 2011) is divided in 4 steps: Acquisition -Reconstruction -Matching -Solving.It shows a complete process of registration with the use of these primitives to reach the required accuracy under specific conditions.One of these conditions concerns the exactness of primitives matching, directly related to the constraints used to solve the registration problem.
In (Rabbani et al., 2007), the search for corresponding objects is done by filtering the possible correspondences using constraints and thresholds.In such an approach, taking the uncertainties of initial approximate transformations into account is not straightforward.Moreover, no information on the probability of a false positive in the established correspondences is delivered; this is the key problem that can ruin the whole process.
Our contributions in this paper are the following: * Corresponding author: Jean-Franc ¸ois Hullo • the definition of an uncertain 4-parameter transformation of 3D lines, and their representation, corresponding to levelled stations of TLS, cf.Section 3; • the description of a probabilistic matching test, using propagation of probability density function resulting of an uncertain transformation, cf.Section 4; • a test based on a real dataset of an industrial scene: its description, results and a discussion, cf.Section 5.

A STATE OF THE ART
TLS data provides 3D features (points, lines, planes) through automatic primitive reconstruction (spheres, cylinders, patches) that can be used for registration: (Rabbani and van den Heuvel, 2005), (Hullo et al., 2011).The matching of 2D features for registration is an active field of research in photogrammetry and computer vision, known as pose estimation problem.Algorithms based on feature points correspondences might be used with precaution in the case of non Euclidean features, such as frames, lines or planes (Pennec, 1999).This assumption is all the more justified when dealing with uncertain primitives, and uncertain transformations of the primitives (Förstner, 2010), (Pennec, 1999).
Complete, unique and minimal representation of 3D lines does not exist.The advances in projective geometry led to a convenient representation of lines, including lines at infinity, the Plücker coordinates.(Förstner, 2010) presents an approach for dealing with uncertain 3D lines and rigid 6-parameters transformations.(Meidow et al., 2009) details the importance of propagation of uncertainties of 2D lines through the entire reasoning chain, from extraction to decision.The application to 3D lines is though not straightforward.
Most laser scanners used for the acquisition of industrial scenes have a bi-axial compensator and can be levelled as well as the  (Roberts, 1988), (Schenk, 2004) and (Bartoli and Sturm, 2005) inspired the representation proposed in 3.2.
For registration, weighted least square resolution can be carried out (Rabbani et al., 2007).Another approach uses probability for solving pose (Chaperon et al., 2011).This approach better taking into account a priori information (uncertainty and coherence).
For matching too, probabilistic frameworks have been developed in the past few years, for example (Tal and Spetsakis, 2010) and (Van Wyk et al., 2004).They pave the way for a better integration of approximated geolocation systems and their uncertainties.

3D LINES, TRANSFORMATIONS AND UNCERTAINTY
3.1 Levelled TLS and 4-parameters transformations We introduce the stations Sti with associated Cartesian coordinate system Ri such as i ∈ I.The origin Oi corresponds to the virtual centre of the laser scanner.R0 is the global, or world, coordinate system.
Considering a levelled TLS, the only directional degree of freedom of the Euclidean transformation τij = τR i →R j from frame Ri to frame Rj is the rotation around the z-axis, with value κij.tij is the 3-vector with coordinates [txtytz]ij representing the translation bringing Oi to Oj.For the frame transformation from Ri to R0, we use the following notation: κi0 = κi and ti0 = ti.Using this formulation, the frame coordinate transformation of a point P from Ri to R0 is : P0 = Mκ i Pi + ti where Mκ i is the rotation matrix around z-axis with value κi.
We use different probability distributions to quantify the uncertainty of these 4 parameters tx, ty, tz and κ.Let t = [txtytz] be a random 3-vector (the translation vector).We introduce a 3dimensional multivariate normal distribution.It can be defined using two parameters (µ, Σ), respectively the mean 3-vector and the 3 × 3 covariance matrix.The probability density function is given by: where |Σ| is the determinant of Σ.
Let κ, the orientation of a station, be a random variable.We introduce a von Mises distribution, or circular normal distribution.It can be defined using two parameters (χ, λ), respectively describing the location and the concentration of the random variable κ.The value λ parameter (λ ∈ [0; ∞]) is 0 when the distribution is uniform, and increases as the distribution concentrates about the angle χ.The probability density function of κ is given in (Mardia, 1972) by: where I0 is the modified Bessel function of order 0. For an order j, this function is given, using the Gamma function (Abramowitz and Stegun, 1964), by: 4) We can now introduce A(λ) that allows us to express the usual linear variance σ 2 of κ: For indoor scenes, several geolocation systems can approximate the 4 parameters κ, tx, ty and tz, such as INS or odometer.The parameters of the a priori distribution of the approximated values are known through hardware specification or calibration of the geolocation system.The uncertainty of the approximated positions and orientation of the stations of our dataset are given in 5.1.

3D lines: representation and transformation
Several representations exist for 3D lines.Some of them are complete (2 points), unique (Plücker) or minimal (two planes intersection, Denavit-Hartenberg).In the case of incomplete representations, we use several representations, called charts, that overlap in order to handle all existing lines.We call atlas a set of charts.
The use of one of these representations depends on the computations and analysis carried out on the lines.In our case, we are interested in the uncertainty of the line coordinates through an uncertain transformation.Statistical analysis requires a minimal representation, as complete as possible.Since we are not dealing with projective geometry, we do not need to take into account lines at infinity.We thus introduce a representation inspired by (Roberts, 1988), illustrated in Figure 2. Frame transformation relations τij = τR i →R j of the coordinates of a line are also presented.The following parameters are used in our representation: 2 is the angle between the z-axis of the frame and the direction vector of the line One can note that almost vertical lines (θ ≈ 0) are a special case in this representation since ϕ is note defined.Since many lines are vertical in industrial scenes, we need to complete the atlas using the following chart for quasi-vertical lines.θ and ϕ are identical to horizontal chart.For defining the position of a quasivertical line, we use (x, y) as the intersection of the line and the O X Y plane of the frame.
For the horizontal lines, the sign of u0 is undefined since ϕ is only determined modulo π.Practically, for matching tests, this case can easily be handled.
There are two main advantages to use this atlas: 1. influences of the rotation and translation transformation parameters are isolated 2. the "singularities" reflect in some way the information that the type of lines delivers for registration: • Quasi-Vertical lines: ϕ is ill-determined.A matched pair of quasi-vertical lines do not provide confident information on the κ or the tz • Oblique line: a pair of oblique lines provides information for both orientation and position parameters, • Quasi-Horizontal line: a pair of horizontal lines provides confident information on κ, (±π/2) and on tz, though it exists an infinity of tx and ty solutions.

Propagation of lines through uncertain approximated transformations
Our dataset is composed of a set Ei of reconstructed lines L i α for each station Sti and approximations of the transformations parameters from each station coordinate system Ri to the global coordinate system R0.The coordinates of lines L i α are considered as certain in the frame Ri of their original station Sti, cf. (5.1).
To be able to compare a line L 1 A reconstructed in St1 to another line L 2 B reconstructed in St2, we have to express them into a common coordinate system.Let assume that we express the coordinates of L 1 A = (θA1, ϕA1, uA1, vA1) into R2, with L 2 A = (θA2, ϕA2, uA2, vA2), using transformation τ12 from R1 to R2.This transformation is composed of two approximated transformations τ10 and τ02 = τ −1 20 .Then L2 A = τ12(L . Since τ10 and τ20 are uncertain and follow a priori distributions, the resulting line L2 A will have uncertain coordinates, following a posteriori distribution.The parameters of the a posteriori distributions of the lines coordinates resulting of an uncertain transformation are presented below.Notations of paragraph 3.1 are used for the distributions parameters.Tilt θA2: as described in (3.2), θ is invariant to the 4-parameters transformation.Then θA1 = θA2.
Direction φA2 of non vertical lines: follows a von Mises distribution with parameters (χ2, λ2).For the mean direction, we have: χ2 = ϕA2 = ϕA1 + κ12 = ϕA1 + κ1 − κ2.The λ2 parameter is given by the equation: where A is the ratio of the modified Bessel functions, (cf.3.1).λκ 2 is estimated using an approximation A inv of the inverse of the A function, given in (Fisher, 1993).Consequently: Position (uA2, vA1) or (xA2, yA2): the distribution of the resulting 2D-random variable corresponds to a Gaussian random variable rotated by a von Mises random variable.For both charts we use as approximation a bi-variate normal distribution with parameters (µ2, Σ2).This approximation is assumed as valid for small variances.The algebraic computation of these parameters can be done developing the two first moments of the resulting random variable.However, in this first implementation, we currently used a Monte-Carlo estimation of these parameters.Probability density function of the bi-variate gaussian is given by: where x = (x, y) ∈ R 2 .

PROBABILISTIC MATCHING OF 3D LINES WITH APPROXIMATED STATIONS 4.1 A priori probability density functions as matching score
A possible approach for computing a matching score between the lines L 1 A and L 2 B is to compute the density of probability of L 2 B using the a posteriori distributions of L 2 A .We can reformulate this idea by the following question:"Knowing the uncertainty of the transformation τ12, what is the probability for the line L 2 B reconstructed in St2 to be identical to L 1 A reconstructed in St1?".We can also express this idea by using the following formulation: , where d is the probability density value of the transformed coordinates of LA.
One of the advantage of the introduced representation is to separate the influences of the orientation and translation parameters of the transformation τ12.Using the formulas and the notations for lines, detailed in 3, we express the probability score for both direction and position of the lines as follows: -Non-vertical lines: -Vertical lines: where g and f are respectively the von Mises and bi-variate normal a posteriori probability density functions of L 1 A .The higher these values are, the more probable is the match between LA and LB.

Steps of the probabilistic 3D lines matching test
Using the previously described uncertainty propagation and density values, the authors proposed a matching test for two 3D lines.They might represent the axes of reconstructed cylinders from two distinct stations of levelled TLS.The steps of this test are illustrated in Algorithm 1.
Two preliminary tests concerns the invariants: radius r and tilt θ.The radius of the cylinders are invariant to any rigid transformation of a cylinder.The tilt is invariant to the 4-parameters transformation defined above.TLS point clouds are currently dense and reconstruction algorithms improved.We can reasonably considered the radius as certain enough to use a threshold, depending on the algorithm and the data set used.If more information is provided on the reconstruction step, one can use a threshold function.
The orientation test is performed on ϕB2 using the a posteriori probability density function of ϕA1.The return s1 is a density value.The second test evaluates the position of the two lines.For the non-vertical ones, the two tested lines are first made parallel, to compare values of (u, v).This comparison is carried out using the probability density functions mentioned above.The return s2 is a density value.s is the probabilistic matching score of the two One can note that any matching algorithm can produce: (a) correct matches, (b) false positive (wrong match) and (c) false negative (missed match).Some registration algorithms will better deal with wrong matches than other but all need as correct matches as possible.Error ratios of our implementation tested on a real dataset will be used to assess the presented approach.

Application to a multi-station scene
Let a scene be composed of i ∈ I stations Sti with αi reconstructed cylinders for each station.In each station, the cylinders are unique, such as any cylinder Ai reconstructed in Sti can only be matched with a unique cylinder Bj reconstructed in Stj.This leads to choose a best matching candidate BJ for Ai, using probabilistic score computed, where sJ = max(sj).A threshold based on confidence interval of the distributions reject weak maxima.
A first remark may be made regarding the multi-station computation of a probabilistic matching score.For any line, the match can: not exist in the dataset (cylinder only reconstructed in one station), be simple (cylinder reconstructed in two stations) or multiple (cylinder reconstructed in more than two stations).it can be assumed that the more stations the cylinder is seen from, the more confident the match should be.This confidence should impact the value of the matching score.Probabilistic scores give access Bayesian theorem and Shannon-Gibbs entropy maximisation, a powerful tool for decision making, cf.(Pennec, 1999) and (Mardia, 1972).

The tested scene
To assess our probabilistic algorithm, we used a real dataset of industrial installations.The scene is a (11 m × 5 m × 3 m ) room. 4 TLS stations were needed to get enough coverage for 3D reconstruction.This scene was chosen because of the large relative errors in the lines position caused by approximated position and orientations of the stations.The other frequently configurations in industrial installations are corridors and large halls, (Hullo et al., 2011).Each per-station point cloud is composed of 25 millions of points.Cylinders have been reconstructed using a region growing algorithm.This algorithm returns the following fitting errors: mean error < 0.002 m , maximum absolute error < 0.01 m.In every stations, we reconstructed 9 up to 16 cylinders, for a total of 53 cylinders (cf. Figure 4)."Real" positions and orientations of the TLS stations have been acquired with a total station.The approximations of the positions and orientations of the stations in a global coordinate system have been randomly generated using known expected errors of current indoor geolocation system (Hullo et al., 2011).Since industrial installations have horizontal grounds and that the instrument heights are measurable, the approximation of the tz will not overcome 10 cm.For tx and ty values, geolocation system provides maximum errors of 15 cm to 50 cm in our study case.The maximum expected error of orientation κ ranges from 2 • to 10 • .Given these error intervals, we carried out 3 tests using parameters of the a priori distributions of the approximated transformations, presented in Table 1.Table 1: A priori distribution parameters of the estimated 4parameters transformation used.σ is the square root of the linear variance.
As mentioned in 3.3, our implementation currently uses Monte Carlo samples to estimate the a posteriori bi-variate normal distributions.We also used the Monte Carlo method to confirm that the bi-variate normal approximation of the resulting random variable was valid.Once matching test is performed on every possible match, our implementation searches for a maximal matching score for each line.In our dataset, we have 1031 score to compute for a maximum of 65 theoretically possible matches.
All the tests have been implemented using R, (R Development Core Team, 2011).

Results
Reference matches have been created by using "real" positions and orientations of the stations.34 matches exists in the scene, 27 non vertical and 7 vertical.A visual inspection of each one as also been carried out.Results of our algorithm are given in Table 2.These first results shows that our implementation, with our dataset and confidence intervals, correctly detect most of the matches.However, vertical lines are not as correctly handled as horizontal ones.More tests using other dataset will be performed to evaluate this trend.As expected, small variances return more correct matches than larger.

Discussion
The existence of false positive and false negative is mainly due to the decision step, currently implemented as a simple maximal score detection.With our dataset, even a small change of the position and orientation of the stations heavily impact the relative distance between lines.In this situation, to properly handle the matching task, one need to take into account the whole structure of the scenes using Bayesian theorem.As mentioned in 4.3, the probabilistic score presented in this article offers such possible process.However, the detected matches and corresponding probabilistic scores computed using the current algorithm could be used as the input for a refinement stage.
The representation of the lines, presented in 2, is close enough to the representation of 4-parameters transformation to allow the distinct analysis of an uncertain one.Regarding the handling of vertical lines, the representation used shows the singularities of these lines when using 4-parameters transformations.The more they approach verticality, the more sensitive they are to a reconstruction error.The Fisher-Snedecor distribution conveys the influence of the value of tilt θ on the determination of ϕ and could be used for quasi-vertical lines, instead of the von Mises distribution.
Otherwise, one of the advantage of the presented probabilistic score is to take into account any a priori knowledge of the uncertainty of the transformation parameters.Then, this score might not be use only with a priori uncertainty of geolocation systems, but can be used during an iterative registration process to check the existing matches.

Considered improvements
Regarding the current implementation, some improvements have been mentioned: • the algebra computation of two first moments of the resulting random variable, presented in 3 is planned, • the use of a F-distribution for quasi-vertical lines might help to handle the direction sensitivity when θ is close to 0, • probabilistic description of the contextual information , for example the probability of intervisibilities using a priori CAD models and integration in a global score, • Bayesian approach, with a priori distributions, for better integration the global structure of the scene by maximising the coherency of the network.

CONCLUSIONS AND FURTHER WORK
In this article, we have presented a probabilistic matching test using propagation of a priori distributions.This test has been applied to 4-parameters uncertain transformation of 3D lines.This configuration corresponds to the matching of the axes of cylinders reconstructed from levelled TLS data, with a priori distribution of transformation parameters provided by indoor geolocation systems.
We introduced a representation of the 3D lines allowing statistical analysis of the influence of an uncertain transformation.Directional and Euclidean statistics were used to model the uncertainty of the transformation parameters and line coordinates.An algorithm using this matching test as been carried out on a real industrial dataset.Experimental parameters of the a priori distributions were chosen to reflect levels of uncertainty of current geolocation systems.First results indicates the validity of the method.
Algebra computations will be carried out to complete the stochastic model.The use of a specific distribution, known as Fisher-Snedecor, is one of the lead considered by the authors is to handle the tilt of the lines and its influence.Information theory and probabilistic graphs are studied for the implementation of a structureawareness algorithm.
Once the presented probabilistic approach is completed, we plan to compare it to other approaches used to match primitives, such as (Rabbani and van den Heuvel, 2005).
Further work is to integrate this matching algorithm in a probabilistic registration solver.The idea is to use the probabilities returned by the matching step to guide the registration and viceversa.Longer-term objective is to shorten on-site time of acquisition by keeping at the same time the required accuracy of 1-2 cm of the final point cloud.

Figure 2 :
Figure 2: A representation for 3D lines.M is the closest point of the line L to the origin O of the frame.

Figure 3 :
Figure 3: Probabilistic scores for matching computed with a posteriori distributions of lines coordinates.

Figure 4 :
Figure 4: Dataset used in our tests.a Cylinders reconstructed from 1 over 4 stations b Axes of cylinders in a common frame using approximate transformations c Ground-truth transformations (total station traverse) are used as reference.
Comparison of the results for the 3 tests performed and the ground truth.For each number, vertical (V) and non vertical (Hz) lines are specified.False negative values indicates missed matches.False positive values indicates wrongly detected matches.
ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume I-3, 2012 XXII ISPRS Congress, 25 August -01 September 2012, Melbourne, Australia total stations.Since a 3D line only has 4 degrees of freedom, we can search for a proper representation to deal with uncertain transformations.Minimal representations used in