Open AccessArticle

Survey of Point Cloud Registration Methods and New Statistical Approach

Jaroslav Marek

^1,*,†

and

Pavel Chmelař

^2,†

Department of Mathematics and Physics, Faculty of Electrical Engineering and Informatics, University of Pardubice, Studentská 95, 532 10 Pardubice, Czech Republic

Department of Electrical Engineering, Faculty of Electrical Engineering and Informatics, University of Pardubice, Studentská 95, 532 10 Pardubice, Czech Republic

Author to whom correspondence should be addressed.

^†

These authors contributed equally to this work.

Mathematics 2023, 11(16), 3564; https://doi.org/10.3390/math11163564

Submission received: 30 June 2023 / Revised: 4 August 2023 / Accepted: 11 August 2023 / Published: 17 August 2023

Download

Browse Figures

Versions Notes

Abstract

The use of a 3D range scanning device for autonomous object description or unknown environment mapping leads to the necessity of improving computer methods based on identical point pairs from different point clouds (so-called registration problem). The registration problem and three-dimensional transformation of coordinates still require further research. The paper attempts to guide the reader through the vast field of existing registration methods so that he can choose the appropriate approach for his particular problem. Furthermore, the article contains a regression method that enables the estimation of the covariance matrix of the transformation parameters and the calculation of the uncertainty of the estimated points. This makes it possible to extend existing registration methods with uncertainty estimates and to improve knowledge about the performed registration. The paper’s primary purpose is to present a survey of known methods and basic estimation theory concepts for the point cloud registration problem. The focus will be on the guiding principles of the estimation theory: ICP algorithm; Normal Distribution Transform; Feature-based registration; Iterative dual correspondences; Probabilistic iterative correspondence method; Point-based registration; Quadratic patches; Likelihood-field matching; Conditional random fields; Branch-and-bound registration; PointReg. The secondary purpose of this article is to show an innovative statistical model for this transformation problem. The new theory needs known covariance matrices of identical point coordinates. An unknown rotation matrix and shift vector have been estimated using a nonlinear regression model with nonlinear constraints. The paper ends with a relevant numerical example.

Keywords:

3D range scanning; automatic space mapping; registration problem; estimates of transformation parameters; ICP algorithm; normal distribution transform; feature-based registration; iterative dual correspondences; probabilistic iterative correspondence method; point-based registration

MSC:

65D18; 62J05

1. Introduction

In the registration problem, data describing a certain object in different coordinate systems is available. This arises quite naturally if we measure points on the object from different positions. We thus obtain so-called clouds of points in individual scans. When solving the registration problem, we search for matching (identical) points in different scans and estimate the transformation parameters. The criterion for estimating unknown parameters is often the minimization of the alignment error. Pioneering methods for solving the registration problem were developed in the 1980s in the field of 2D image matching. This important problem of computer graphics crystallized, especially during the processing of medical images. Research gradually moved towards the application of 3D point cloud registration. Applications in object reconstruction, non-contact inspections, medical and surgery support, robot navigation, and autonomous vehicle navigation gradually appeared. An overview of methods used in 2D registration is published in the book [1]. The paper [2] contains a literature survey on the use of the registration problem in robotics and robot navigation. The state-of-the art of registration problem with an emphasis on radar data is dealt with in the book [3].

Nowadays, laser scanning is finding huge use. A comprehensive review of terrestrial laser scanner (TLS) point cloud registration methods is given in ref. [4]. The book [5] expressed the TLS as follows: “TLS systems enable documentation of the whole measured object with all constructional elements. The high scanning speed of the current scanners enables a significant reduction in the time necessary for measurement and, alternatively, an increasing amount of information about the measured object. More data (several million points) enable the creation of complex models or documentation of objects in the form of a detailed point cloud”. A relatively large potential for the use of such models has arisen. As such, laser scanning finds application in various fields: industry, robotics, process control, quality control, architecture, civil engineering, archeology, geosciences, medical sciences, and dentistry. Autonomous laser scanning systems also offer the prospect of saving human lives. They are deployed in dangerous environments such as unexplored areas, abandoned places, or mines. Such systems offer the possibility of getting detailed maps of the unknown surrounding environment without human intervention. In these systems, solving registration problems is crucial. Cf. [5].

Our measurement was made by using an optical rangefinder, which uses an optical filter to spread the green laser beam in a vertical line. See Figure 1. At each detected point of this line, it is possible to determine a 3D point in the scanned environment. Figure 2 demonstrates laser scanning of an indoor room. The optical rangefinder’s measurement head is rotating in 360 degrees around the Z axis to scan the whole surrounding area. The resulting point clouds form two positions and are presented in Figure 3. Using a camera, you can create a colored point cloud.

In this paper, we will mainly discuss methods applicable to 3D registration. In doing so, we will focus on the registration of building objects. Firstly, we present an overview of the most commonly used methods. In addition, we will propose a statistical algorithm that can become an appendix to standard algorithms. In the model, we will assume that the transformation is formed only by shift and rotation in 3D space. Such a transformation is sometimes referred to as a rigid transformation. In this new algorithm, we will consider measurement errors in individual scans. This regression model will allow the synchronous estimation of unknown transformation parameters from all scans. The prerequisite for using this approach will be the knowledge of identical points. These points can be earned through standard registration methods. In special geodetic problems, points on the object can be given by placing reflective targets, which are then focused in individual scans. Here, it is not necessary to look for identical points, as they are determined by reflection points.

Let us return to the description of an example of a typical 3D registration problem. Measured coordinates of points are given in the point clouds, which can be obtained from different scans from different positions. Figure 3 gives a point cloud registration example. The registration problem’s aim is to create a compact point cloud with all possible details.

The structure of the article is as follows. In the introductory chapter, the framework of existing knowledge was outlined, the assumptions and goals of the proposed statistical model were described, and the realization of a specific laser measurement was presented. The second section seeks to examine and discuss the main currents of thought in the theoretical and practical solution of registration problems without attempting to present the details of all the methods. In the third chapter, we will state a new statistical method, which exploits a regression model with constraints. In the fourth chapter, a concluding discussion is made.

A fundamental shortcoming of commonly used methods is the fact that the covariance matrices of the measured points are often not used or are replaced by identical matrices. Further, the estimate of the variance matrix of the estimated coordinates is not known. We find this estimate using an indirect measurement regression model with a system of nonlinear constraints. Separately, this model is applicable in the case of geodetic measurements, when identical points are marked with reflective targets. After gaining pairs of identical points from commonly used methods, the new approach will find application in the process of object modeling.

1.1. Optical Rangefinder

The principles of the 3D technique for object description are based on the triangulation method. The idea of measuring the distance and then determining the point cloud is mentioned in the article [6]. A rotating tripod carried the entire measuring device; see Figure 1. A Basler camera with an output of 2590 × 1942 pixels was used to capture the image in high quality. The output power of 200 mW of the used laser diode made it possible for good recognition after vertical swapping. An optical filter with a 90-degree angle scattered the laser spot. A powerful and fast-stepping motor enabled the measuring head to rotate 360 degrees. Captured colored images are processed by image processing methods to get precise laser line determination in the individual scanned images. For each detected point on the line, the 3D position in the surrounding scanning area is calculated. Image information is also used in the point cloud creation process to get a colored point cloud; see ref. [6].

1.2. 3D Point Cloud

The Figure 2 left side shows an example of the implementation of one frame measurement. Laser line segmentation was performed using a green-colored Gaussian Mixture Model (GMM). For more details, see the article [6]. By using laser mask analysis, 10 different distances were measured. The total number of N measured points in one image was 80. The mobile device is located at the coordinates

P = (P_{X}, P_{Y}, P_{Z})

. The point P is the origin of one automatic 360-degree scan. Through measurement data processing, the 3D point cloud is created. An example of a 3D point cloud is in Figure 3. The used coordinate system is Cartesian, and each point in this point cloud is described by a vector

p = (p_{X}, p_{Y}, p_{Z})

, where

p_{X}

p_{Y}

, and

p_{Z}

are coordinates in the appropriate axis.

The black dot marks mobile platform position

P_{1}

. This point cloud covers the measurement from one position, and there are some important points missing in the top corner. The whole scanning process’s intention is to create a compact map of a scanning space. For this purpose, it is necessary to provide the measurement from more positions in a scanning space. For changing the measurement position, the used measurement device can be, for example, placed manually in a new position or mounted on a mobile device, which provides precise information about its position in space. We place a mobile platform manually in a scanning environment. The differences in position we measure with the precise laser rangefinder.

2. A Survey of Methods

In the following subsections, we will present the ideas of the most commonly used methods developed to solve the registration problem.

2.1. ICP Algorithm

ICP is the oldest registration algorithm; cf. [7,8]. If the user’s measurement or device’s localization is accurate, the registration can be determined directly by separate measurements. Figure 3 describes the connection to one group. However, imprecise measurement sensors and erroneous self-localization have to be assumed for registration. The frequently used method is called Iterative Closest Points (ICP); see refs. [9,10].

Let’s have two independently gained sets of 3D points, the measurement from position

P_{1}

and

P_{2}

, we want to estimate the transformation

(R, t)

forming by a rotation matrix

R

and a shift vector

t

which minimizes the cost functional

\begin{matrix} E (R, t) & = & \sum_{i = 1}^{N_{m}} \sum_{j = 1}^{N_{d}} w_{(} i, j) | | {\hat{m}}_{i} - (R {\hat{d}}_{j}) + t {| |}^{2} . \end{matrix}

(1)

The parameter

w_{i, j}

takes the value 1 when the i-th point of

\hat{M}

determines the same 3D position of the point as the j-th point of

\hat{D}

. If it specifies another point, then

w_{i, j}

is 0. Next, two parameters need to be determined: the corresponding points and the function

E (R, t)

that minimizes

E (R, t)

with the use of corresponding points. The ICP algorithm iteratively computes the point correspondences. In each step of the algorithm, the algorithm searches the closest points as correspondences and computes the transformation

(R, t)

for the minimizing value

E (R, t)

2.1.1. ICP Algorithm Results

The ICP algorithm application result is shown in Figure 4, and the final RMS value between both point clouds is 0.167. The point cloud contains points from the second measurement that fill the empty area in the first measurement (top corner). The new point cloud contains more measurement points and gives more details about the scanned room. In Figure 4, on the left side are both registered point clouds from Figure 3. On the right side, there is a colored version.

2.1.2. ICP Algorithm Problem

By using the ICP method, there are two main problems. The ICP algorithm is point-based, and it does not take into account the local surface and shape around each point. The ICP algorithm is also time-consuming in searching for the nearest neighbor. Even if an efficient search data structure such as a kd-tree is used, the estimation time for a large number of points is long. If both scans are with a bigger distance offset and rotated in all axes by about several tens of degrees, the registration result is obtained with error.

2.2. Normal Distribution Transform (NDT)

The method was proposed in the publication [11]. The basic pillar of this algorithm is the capture of data by a statistical model. The method does not focus on the simple use of individual measurement points but instead constructs a model of mixtures of normal distributions and examines the probabilities of the surface point occurring in a certain position; see ref. [12]. This approach provides a piece-wise smooth introduction to the point cloud model with continuous first- and second-order derivatives. This representation is followed by standard numerical optimization methods for the registration problem. The NDT algorithm consists of several steps: firstly, the whole point cloud is dismembered into small subdomains (cells), and for each cell, statistical characteristics such as the mean and the covariance matrix are computed. We consider that the coordinates of the reference scan surface points have a D-dimensional normal distribution. The likelihood function is

\begin{matrix} p (x) = \frac{1}{{(2 π)}^{\frac{D}{2}} \sqrt{| Σ |}} exp (- \frac{1}{2} {(x - μ)}^{'} Σ (x - μ)), \end{matrix}

(2)

where

μ

is the vector of mean values and

Σ

is the covariance matrix of the reference points on the scanned surface in the cell in which

x

is located. The constant

c_{0}

denotes the expression

{(2 π)}^{\frac{D}{2}} \sqrt{| Σ |}

. The estimates of the mean and covariance matrix are

\begin{matrix} \bar{μ} & = & \frac{1}{m} \sum_{k = 1}^{m} y_{k}, \end{matrix}

(3)

\begin{matrix} \bar{Σ} & = & \frac{1}{m - 1} \sum_{k = 1}^{m} (y_{k} - μ) {(y_{k} - μ)}^{'}, \end{matrix}

(4)

where

y_{k}

k = 1, \dots, m

are coordinates of reference scan points form a cell. According to the eigenvalues of the covariance matrix, it is possible to estimate the shape of the cell. This will help in the point cloud registration process. In the NDT method, we look for such a position of the actual scan that when the maximum probability occurs, the points of the actual scan are located on the surface of the reference scan. Searched Helmert transformation parameters (i.e., translation and rotation) offer a pose estimate in the given scan; for encoding, we will use the vector variable p. In the current scan, we have point clouds

X = {x_{1}, \dots, x_{n}}

The displacement of the point

x

in space by the position

p

is described by the spatial transformation function

T (p, x)

. We consider the best estimate of the true position

p

to be the one that maximizes the likelihood

\begin{matrix} L & = & \prod_{k = 1}^{n} p (T (p, x_{k})) . \end{matrix}

(5)

For 3D-NDT registration, the Euler angles are useful because of their natural advantages. This representation does not need a numerical optimization procedure with easier-to-calculate derivatives. If we use the representation using Euler angles, there are six transformation parameters to estimate: three unknown translation parameters and three unknown rotation parameters. The rotation matrix can be expressed by Cardan angles and angles of yaw, pitch, and roll, ref. [13]. The pose can be explained using the six-dimensional vector

p_{6} = {[t_{x}, t_{y}, t_{z}, ϕ_{x}, ϕ_{y}, ϕ_{z}]}^{'}

using the Euler angles. Then the 3D transformation is given by formula:

\begin{matrix} T_{E} (p_{6}, x) = R x + t = R_{x} R_{y} R_{z} x_{k} + t, \\ R = (\begin{matrix} c_{x} c_{z} & - c_{y} s_{z} & s_{y} \\ c_{x} c_{z} + s_{x} s_{y} c_{z} & c_{x} c_{z} - s_{x} s_{y} c_{z} & - s_{x} c_{y} \\ s_{x} s_{z} - c_{x} s_{y} c_{z} & c_{x} s_{y} s_{x} + s_{x} c_{z} & c_{x} c_{y} \end{matrix}), \\ t = {(t_{x}, t_{y}, t_{z})}^{'} \end{matrix}

(6)

where

c_{i} = cos ϕ_{i}

and

s_{i} = sin ϕ_{i}

2.3. Feature Based Registration

The feature-based matching extracts distinguishing features from the range of images and uses corresponding features to calculate scan alignment. For most robots, the detection of closed loops is realized from camera data [14]. The most common functions of the detector are: SURF; SIFT; GLOH; Shape Context; PCA; Moments; Cross-correlation; and Steerable filters [9]. A major aspect of feature detectors is their invariance with individual transformation parameters: scale, shift, rotation, and illumination changes. A comprehensive survey of detectors and descriptors is given in [15]. These methods use the Gaussian function for image spectral analysis

\begin{matrix} G (x, y, σ) & = & \frac{1}{2 π σ^{2}} exp (\frac{x^{2} + y^{2}}{2 σ^{2}}), \end{matrix}

(7)

where x and y are pixel coordinates in an image, and

σ

is the Gaussian parameter. The n-times filtration provides scale invariance. The input images of the convolution with variable scales give the differences of Gaussian (DoG), which are the image features. For a significant increase in computation speed, the SURF method uses Gaussian function approximation. The feature-based navigation system is presented in ref. [16]. The autonomous robot is guided along a given path. The navigation is based on exteroceptive sensors, and a simple, effective mathematical model is suitable.

2.4. IDC

The iterative dual correspondence (IDC) algorithm [17] can be considered an extension of ICP. The main mission is to accelerate the convergence of the rotation matrix when estimating position in scan registration. In this method, two rules are introduced to determine the correspondence. The first is similar to the ICP. In each iteration, the rotation and translation pair

τ_{1}

are determined using the closest points. The second criterion, called ’the matching range-point rule’ deals with the point selection of a new set of corresponding points. In this criterion, corresponding points are identified within the permissible angular interval. An equation in 2D space has the form

[φ - t_{φ}, φ + t_{φ}]

t_{φ}

and is limited (how far the algorithm should seek). The matching range-point metric is determined as follows:

corresponding (x) = {\arg \min}_{x^{'}} (| r - ρ |)

where

x = [φ, r]

and

x^{'} = [ϕ, ρ]

, and

ϕ - t_{ϕ} \leq ρ \leq ϕ + t_{ϕ}

The corresponding point is located within the specified angular limits. In 3D space, an interval is a cuboid-shaped “window” whose projections into the planes are rectangles. A three-dimensional point

x

is given by polar coordinates

[ϕ, θ, r]

, where

ϕ

and

θ

are the latitude and longitude angles. This window is bound from

[ϕ - t_{ϕ}, θ - t_{θ}]

[ϕ + t_{ϕ}, θ + t_{θ}]

. A comparison of the quality of the solution to the registration problem led to the conclusion that the IDC algorithm is more robust than the ICP. However, with a large bias in the initial estimate of the position compared with the actual value, greater inaccuracies in the solution occurred more with the IDC than with the ICP. The IDC method is much more sensitive to a biased initial solution. Cf. [17].

2.5. pIC

Another approach is a registration algorithm called the probabilistic iterative correspondence method (pIC), cf. [18]. In this algorithm, scan noise data is combined with an initial estimate of position. The coordinates of the scan points and the initial estimate of the position are assumed to be random variables with zero mean Gaussian noise. For the best result from matching, the covariance matrix is based on the sensor positions and robot odometry. In the beginning, points from both scans are selected to create a subset composed of statistically compatible points. The construction of the pairing fitness function is based on the minimization of the Mahalanobis distance

d_{m} (x, y)

. of point

x

and

y

\begin{matrix} d_{m}^{2} (x, y) & = & | | x - y | | Σ^{- 1} | | x - y | | \end{matrix}

(8)

where

Σ

is a known estimate of the covariance matrix, including sensor properties and system conditions. If the Mahalanobis distance does not exceed the specified threshold, these points are considered compatible. The set of points satisfying this constraint formulates a subset of A. The method for determining the minimum argument of the fitness function, i.e., estimating the point

y_{i} \in A

that best corresponds to the point

x

, is the integration of each admissible locations

x

and each admissible sensor coordinates considering their Gaussians estimates.

The results of the pIC algorithm were compared with those of IDC and ICP. pIC estimates were more robust with respect to a biased initial solution.

The pIC method shows faster convergence, about 25 %. Cf. [18].

2.6. Point-Based Probabilistic Registration

In this method, the measured coordinates of the discrete points of the reference scan are replaced by probabilistic functions, cf. [19]. The basis for estimating the probability of a scan point from the current cloud is to trace a ray rising from the estimated position, along the direction associated with each measurement from the current scan, to the nearest surface in the reference scan. Triangulation of the reference coordinate system scan is then performed to create an approximation of the surface. The estimated range of this measurement is determined by the length of the beam.

The likelihood function of a scan point is determined from a mixture of Gaussian values that are centered in an estimated range, with variance determined by known scanner characteristics, and a uniform distribution that is also determined with respect to scanner uncertainties. The probability of the distance r obtained by measurement (the distance to the cloud point

x

) depends on the theoretical distance (the mean of the Gaussian value). This value is calculated through the mixture model using the Euclidean metric between scan points

x

and the triangulated reference scan surface. The calculation uses the expected metric

d_{e} (x, Y, p)

pursuant to the direction of the beam

x

, the reference cloud Y, and the position

p

. The probability that the position of the current cloud is

p

is given by multiplying probabilities

\prod_{x \in X} p (x | d_{e} (x, Y, p))

The algorithm tries to optimize the likelihood function. The presented results by the authors include a pair of building scans in different coordinate systems, and it is discussed that the method gives results with less uncertainty in comparison with the ICP on this data set. Cf. [19].

2.7. Gaussian Fields

The next developed registration method has the name Gaussian fields, cf. [20]. This approach is similar to the Normal Distribution Transformation (NDT). A Gaussian mixture model (GMM) becomes the main tool here for measuring the scan’s statistical characteristics, the spatial distance between the points, and the surface similarity around the points. The comparison of points is realized in a multidimensional space, the dimension of which is given by the spatial dimension of the coordinate system, which is increased by the number of attributes. The proximity and similarity of points

x

and

y

from different clouds is

\begin{matrix} F (x, y) & = & exp (- \frac{| | x - y | |}{d_{g}^{2}} - {(S (x) - S (y))}^{'} \times \\ D^{- 1} (S (x) - S (y))), \end{matrix}

(9)

where

S (x)

is the 3D-moment shape characteristics of the surface passing through

x

. Equation F(x, y) describes a Gaussian function centered at

y

. The parameter

d_{g}

describes the decomposition rate according to the spatial function, and the diagonal matrix

D

determines the penalty given by differences in attributes. The fitness function is

\sum_{i, j} F (x_{i}, y_{i})

. Cf. [21].

2.8. Quadratic Patches

In this method, the reference scanned surface is expressed implicitly when the approximant is constructed using the square of the distance from the surface. The minimization of the distance metric given by the sum of squared distances defines the criterion to connect points from the current scan. The local shape of the surface is fitted by approximants using quadratic terms, i.e., the second-order derivative. The optimum of this surface representation can be found by Newton’s optimization method. Cf. [22].

2.9. Likelihood-Field Matching

This variant of registration works with similar principles as NDT. The method is often used by developers for 2D sonar scan registration. The large uncertainty of the angular resolution of sonars, compared with lidars, brings complications to the successful identification of registration points. In the method, a sequence of scans is first accumulated using robot odometry to obtain multi-point scans. Such an additional procedure can be applied to the input data for all scan registration algorithms. Thus, a new group of methods was created with the following prefix: sNDT, sICP. The method, named LF/SoG (short for probability field introduced as a sum of Gaussians), is very comparable to NDT when the points in one scan are corrected to a normal distribution according to the points from the next scan. In this method, an identity matrix is used instead of the covariance matrix estimates. Instead of a discrete explicit grid, the entire reference scan is passed through the circular window by replacing the values of the internal points in the sliding window with their centers of gravity. A Gaussian approximant is centered at the point of the resampled cloud, with the mean given by the point pose and the covariance matrix set to the identical matrix. In the process of registration, not only the nearest points are used for the calculation of the optimization criterion but also points at a defined distance from the Gaussian function. Cf. [23].

2.10. CRF Matching

Conditional Random Fields (CRF) is a probabilistic approach used in the construction of relational information models. Unlike hidden Markov models or Markov random fields, the use of CRF does not require independence of observations. Originally, CRFs were implemented in the field of computational linguistics, where they were used to label sentences according to parts of speech. Conditional random fields have been implemented to register 2D scans. In the algorithm, the CRF created for each point in the current scan introduces a hidden node. In turn, each hidden node is assigned to its corresponding point in the reference scan. The hidden node is also linked to a data node that contains certain functions defined for each scan point. On the selected training data, the model learns and searches for suitable parameters. Cf. [24].

2.11. Branch-and-Bound Registration

A branch-and-bound strategy is another registration algorithm presented in [25]. First, the translational part of the position space is discretized in several resolutions with ordered levels from the coarsest to the finest. The sequence of poses is analyzed at a certain level of the hierarchy. At the second lowest hierarchical level (branching), the value of the score is determined, on the basis of which sub-nodes are eliminated (boundary). This application is often used in forestry. The application is suitable for highly unstructured environments. It is particularly interesting in 2D analysis. The branch-and-bound procedure is used in the translation part of the transformation. This method is used in applications aimed at modeling the orientation of the robot based on various sensors. The initial value of the rotation matrix plays an important role here. In a space with six degrees of freedom (position in a space and orientation angles yaw, pitch, and roll angles), the optimization is complicated. Cf. [25].

2.12. Registration Using Local Geometric Features

Compared with the point algorithms described earlier, registration can also be implemented using more descriptive local geometric elements [26]. An ideal local property descriptor is invariant to rigid motion so that the corresponding surface can be found regardless of initial scan positions. The correspondence problem is solved when sufficiently prominent features are found. Furthermore, a global adaptation of the surface is possible, i.e., the realization of registration without an initial estimate of the position. The most commonly used surface description method is ’spin-images’. The spin-image is constructed at an oriented point, which is a surface point with a specified normal vector. Another alternative gives ’surface signatures’. These are similar to spin images but use surface curvature instead of point density. Another alternative surface description is the ’splash’. Spatters are determined around oriented points and involve capturing the point cloud using triangulation methods. All these description techniques are used to express the geometric features of both scans. The main aim is to calculate the estimates of rotation and translation parameters between the scans in such a way that we will obtain their minimal mutual distance. Cf. [26].

2.13. PoitnReg

Another registration approach for a wide range of three-dimensional scans is presented in ref. [27]. This paper deals with the pitfalls of 3D range scans. The rangefinder translates coordinates to a chosen coordinate system, and each scan origin is translated by the X, Y, Z coordinates. The data from an RTK GPS receiver placed directly above the scanner position is available for measurement. The PointReg algorithm evaluates the rotation combinations about the Z-axes of two scans to estimate the rotation pair with the lowest RMS error. The selected points from the first scan are compared with triangles formed by the nearest three points in the second scan.

This algorithm’s initial stage needs interaction by the user to ensure that the registration is determined only at points located on the object. Arbitrary moving objects in the scans are eliminated to ensure low inconsistencies between acquired points. The presented result showed that the algorithm needed only a few integrations, and for the wide-range scan, the PointReg algorithm overcomes the standard ICP.

2.14. The Automatic Mapping of Parametric Objects (Buildings) from Indoor Point Clouds

The recent paper [28] focuses on sophisticated data processing, modeling, and data visualization after the registration process. The point cloud registration process is done by high-quality rangefinder software. This research aims at the volumetric and the parametric building model, which additionally incorporates contextual information such as global wall connectivity. Unlike other methods working with the surface of the monitored object, this reconstruction problem works with the representation of objects/buildings using parametrically described and interconnected volume elements. The automatically reconstructed walls between adjoining rooms from opposite surfaces of vertical walls were extracted in the input clouds, and at the present discretion, the globally satisfactory connectivity of all elements. At the same time, with a precise mapping of wall thickness, the result is a complex volumetric approximation of an object’s wall elements. The described algorithm deals with many challenges, including removing outliers, automatic wall width estimation, calculating individual room area, window and door detection, and exit path determination for the scanned building.

The resultant model is easy to edit in software, and there can be simulated building changes, for example, in the case of future reconstruction. The presented approach was successfully tested on several real building measurements.

2.15. Basic Registration Methods Comparison

In Table 1 and Table 2, the basic method comparison from Section 2.1, Section 2.2, Section 2.3, Section 2.4, Section 2.5, Section 2.6, Section 2.7, Section 2.8, Section 2.9, Section 2.10, Section 2.11, Section 2.12, Section 2.13 and Section 2.14 is summarized. Individual algorithms are compared based on individual points or some space modeling, handling covariance matrix, important requirements, their speed processing, computation complexity, precision and advantages, necessity of human intervention, suitable application, and their year of origin.

Due to various scanning sensors’ errors, different registration methods can fit the concrete algorithm in Table 1 and Table 2. Generally, the statement can be valid, if measurement data are precise and point-based algorithms can perform better. The methods that use some space modeling can be more suitable for noisy input data. The concrete application area can also be important in the algorithm selection because some algorithms are specific to concrete applications. Despite the number of mentioned algorithms in Section 2. the main standard in 3D range scanning is mainly the ICP and NDT and their deviates. They are used and implemented in point cloud processing libraries. The method from Section 2.14 is not mentioned because it is a complex handling of building models, and it is not impossible to compare it with others.

3. The Optimal Statistical Transformation Model Taking into Account Errors in Coordinate Systems

This method can be added as an appendix to any registration algorithm presented in Section 2. The registration algorithm will provide data on the found identical points, from which the transformation parameters and their covariance matrices are subsequently estimated in the regression model.

The statistical model consisting of two scans has been studied in ref. [29], under the assumption that stochastic coordinates are given only at the second coordinate system and the points of the first scan are error-less.

A more complex problem with 3D transformation arises when measurement errors are respected in both coordinate systems. Further, the following notation is used.

Y_{i}

…

3 (n_{i} + n_{i + 1})

-dimensional random vector of realizations, which are obtained by the data in the i-th position of the laser sensor,

i = 1, 2, 3, 4

(in the following text four device positions will be considered),

n_{i}

… number of points P in the i-th laser sensor position; (if

i = 4

, then instead of

i + 1

is the number 1).

The i-th device position must be chosen so that the i-th and

(i + 1)

-th measured object side must be covered by measurement. The often arbitrary i-th position may be selected as the basis for the coordinate transformation into a chosen system of coordinates.

If the chosen coordinate system is equivalent to the first cloud system then the expected value of vector

Y_{1}

\begin{matrix} E (Y_{1}) = (\begin{matrix} {fi}_{1} \\ β_{2} \end{matrix}), β_{1} = {(β_{1, 1}^{'}, \dots, β_{1, n_{1}}^{'})}^{'}, \\ β_{2} = {(β_{2, 1}^{'}, \dots, β_{2, n_{2}}^{'})}^{'}, \end{matrix}

(10)

where

β_{1, j}

are three-dimensional coordinates of the j-th point P on the first side of the object and

β_{2, j}

are three-dimensional coordinates of the j-th point on the other side of the object. If

i \neq 1

, then

\begin{matrix} E (Y_{i}) = (\begin{matrix} α_{i} \\ α_{i + 1} \end{matrix}), \end{matrix}

(11)

where instead of

i + 1

for

i = 4

is the number 1.

3 n_{i}

-dimensional vector

α_{i}

includes the coordinates of

n_{i}

points on the i-th object side at the i-th device coordinate system. Similarly,

α_{i + 1}

has a length

3 n_{i + 1}

when this vector includes the coordinates of

n_{i + 1}

points on

(i + 1)

-th object side at i-th device coordinate system. It can be considered

\begin{matrix} α_{i} & = & 1_{n_{i}} \otimes γ_{i} + (I_{n_{i}, n_{i}} \otimes T_{i}) β_{i}, \end{matrix}

(12)

\begin{matrix} α_{i + 1} & = & 1_{n_{i + 1}} \otimes γ_{i} + (I_{n_{i + 1}, n_{i + 1}} \otimes T_{i}) β_{i + 1}, \end{matrix}

(13)

where

β_{i}

β_{i + 1}

are vectors of point coordinates on i-th side or

(i + 1)

-th object side respectively. The 3D-vector

γ_{i}

and the orthogonal matrix

T_{i}

are unknown and must be appropriately estimated.

Therefore

\begin{matrix} E (Y_{i}) & = & (\begin{matrix} 1_{n_{i}} \otimes γ_{i} \\ 1_{n_{i + 1}} \otimes γ_{i} \end{matrix}) + \\ (I_{n_{i} + n_{i + 1}, n_{i} + n_{i + 1}} \otimes T_{i}) (\begin{matrix} β_{i} \\ β_{i + 1} \end{matrix}) . \end{matrix}

(14)

Further

\begin{matrix} T_{i} = (\begin{matrix} T_{i, 11}, & T_{i, 12}, & T_{i, 13} \\ T_{i, 21}, & T_{i, 22}, & T_{i, 23} \\ T_{i, 31}, & T_{i, 32}, & T_{i, 33} \end{matrix}), \\ T_{i, j} = (\begin{matrix} T_{i, 1 j} \\ T_{i, 2 j} \\ T_{i, 3 j} \end{matrix}), i, j = 1, 2, 3 . \end{matrix}

(15)

For simplicity, the situation will be discussed further on the model with 4 scans.

3.1. Model with Four Scans

The model with four scans is described as follows

E (\begin{matrix} Y_{1} \\ Y_{2} \\ Y_{3} \\ Y_{4} \end{matrix}) \sim [(\begin{matrix} β_{1} \\ β_{2} \\ α_{2}^{(2)} \\ α_{3}^{(2)} \\ α_{3}^{(3)} \\ α_{4}^{(3)} \\ α_{4}^{(4)} \\ α_{1}^{(4)} \end{matrix}), Σ],

(16)

\begin{matrix} Σ & = & (\begin{matrix} Σ_{1}, & 0, & 0, & 0 \\ 0, & Σ_{2}, & 0, & 0 \\ 0, & 0 & Σ_{3}, & 0 \\ 0, & 0 & 0, & Σ_{4} \end{matrix}) \end{matrix}

(17)

and constraints (for

s = 2, 3, 4

) are

\begin{matrix} {T^{'}}_{s, 1} T_{s, 1} = 1, {T^{'}}_{s, 1} T_{s, 2} = 0, {T^{'}}_{s, 1} T_{s, 3} = 0, \\ {T^{'}}_{s, 2} T_{s, 2} = 1, {T^{'}}_{s, 2} T_{s, 3} = 0, {T^{'}}_{s, 3} T_{s, 3} = 1, \end{matrix}

(18)

and

\begin{matrix} α_{2}^{(2)} - 1_{n_{2}} \otimes γ_{2} - (I_{n_{2}, n_{2}} \otimes T_{2}) β_{2} & = & 0, \\ α_{3}^{(2)} - 1_{n_{3}} \otimes γ_{2} - (I_{n_{3}, n_{3}} \otimes T_{2}) β_{3} & = & 0, \\ α_{3}^{(3)} - 1_{n_{3}} \otimes γ_{3} - (I_{n_{3}, n_{3}} \otimes T_{3}) β_{3} & = & 0, \\ α_{4}^{(3)} - 1_{n_{4}} \otimes γ_{3} - (I_{n_{4}, n_{4}} \otimes T_{3}) β_{4} & = & 0, \\ α_{4}^{(4)} - 1_{n_{4}} \otimes γ_{4} - (I_{n_{4}, n_{4}} \otimes T_{4}) β_{4} & = & 0, \\ α_{1}^{(4)} - 1_{n_{1}} \otimes γ_{4} - (I_{n_{1}, n_{1}} \otimes T_{4}) β_{1} & = & 0, \end{matrix}

(19)

(

α_{i + 1}^{(i)}

comprise measured coordinates on the

(i + 1)

-th side of object at the i-th device coordinate system).

Unknown parameters are the following vectors

\begin{matrix} β_{1}, β_{2}, β_{3}, β_{4}, α_{2}^{(2)}, α_{3}^{(2)}, α_{3}^{(3)}, α_{4}^{(3)}, α_{4}^{(4)}, α_{1}^{(4)}, \\ γ_{2}, γ_{3}, γ_{4}, T_{i, 1}, T_{i, 2}, T_{i, 3}, i = 2, 3, 4 . \end{matrix}

The unknown parameter count in the model is

3 (n_{1} + n_{2} + n_{3} + n_{4}) + 6 (n_{3} + n_{4}) + 3 (n_{1} + n_{2}) + 9 + 27

, and the number of constraints is

3 n_{1} + 6 (n_{3} + n_{4}) + 3 n_{2} + 18

. The measurements count plus constraints minus parameters is

3 n - 18

3.2. Linearization of a Model with 4 Scans

Now the linearized model will be described.

The following notation and terminology will be used:

\begin{matrix} Θ_{1} & = & (β_{1}^{'}, β_{2}^{'}, {(α_{2}^{(2)})}^{'}, {(α_{3}^{(2)})}^{'}, {(α_{3}^{(3)})}^{'}, {(α_{4}^{(3)})}^{'}, \\ {{(α_{4}^{(4)})}^{'}, {(α_{1}^{(4)})}^{'})}^{'}, \end{matrix}

(20)

\begin{matrix} Θ_{2} & = & (β_{3}^{'}, β_{4}^{'}, γ_{2}^{'}, γ_{3}^{'}, γ_{4}^{'}, {[vec (T_{2})]}^{'}, \\ {{[vec (T_{3})]}^{'}, {[vec (T_{4})]}^{'})}^{'}, \end{matrix}

(21)

\begin{matrix} Y & = & (Y_{1}^{'}, Y_{2}^{'}, Y_{3}^{'}, Y_{4}^{'}), \end{matrix}

(22)

\begin{matrix} Σ & = & (Σ_{1}, Σ_{2}, Σ_{3}, Σ_{4}) . \end{matrix}

(23)

Let

Θ_{1}^{(0)}, Θ_{2}^{(0)}

are approximate values of parameters

Θ_{1}, Θ_{2}

The linear version of constraints is

b + B_{1} δ Θ_{1} + B_{2} δ Θ_{2} = 0

, where

δ Θ_{1} = Θ_{1} - Θ_{1}^{(0)}

δ Θ_{2} = Θ_{2} - Θ_{2}^{(0)}

and matrices

B_{1}

and

B_{2}

are given in Table 3.

3.3. The Estimates of Unknown Parameters in a Model with 4 Scans

The indirect measurement model with II type constraints is

Y \sim (X Θ_{1}, Σ), b + B_{1} Θ_{1} + B_{2} Θ_{1} = 0 .

If rank

(X_{n, k}) = k_{1} < n,

rank

(B_{1, (q, k_{1})}, B_{2, (q, k_{2})}) = q < k_{1} + k_{2},

rank

(B_{2}) = k_{2} < q,

and

Σ

is a positive definite, then the model is regular.

black of vector

(\begin{matrix} Θ_{1} \\ Θ_{1} \end{matrix})

is given by formula

\begin{matrix} (\begin{matrix} {\hat{\hat{Θ}}}_{1} \\ {\hat{\hat{Θ}}}_{2} \end{matrix}) & = & - (\begin{matrix} {(X^{'} Σ^{- 1} X)}^{- 1} B_{1}^{'} Q_{1, 1} \\ Q_{2, 1} \end{matrix}) b + \\ + (\begin{matrix} I - {(X^{'} Σ^{- 1} X)}^{- 1} B_{1}^{'} Q_{1, 1} B_{1} \\ - Q_{2, 1} B_{1} \end{matrix}) {\hat{Θ}}_{1}, \end{matrix}

(24)

where

{\hat{Θ}}_{1} = {(X^{'} Σ^{- 1} X)}^{- 1} X^{'} Σ^{- 1} Y

(estimate non-respecting the parameters

Θ_{1}, Θ_{2}

constraints); covariance matrix is

\begin{matrix} var (\begin{matrix} {\hat{\hat{Θ}}}_{1} \\ {\hat{\hat{Θ}}}_{2} \end{matrix}) = (\begin{matrix} var ({\hat{\hat{Θ}}}_{1}), & cov ({\hat{\hat{Θ}}}_{1}, {\hat{\hat{Θ}}}_{2}) \\ cov ({\hat{\hat{Θ}}}_{2}, {\hat{\hat{Θ}}}_{1}), & var ({\hat{\hat{Θ}}}_{2}) \end{matrix}), \end{matrix}

(25)

where

\begin{matrix} var ({\hat{\hat{Θ}}}_{1}) & = & {(X^{'} Σ^{- 1} X)}^{- 1} - {(X^{'} Σ^{- 1} X)}^{- 1} B_{1}^{'} \times \\ Q_{1, 1} B_{1} {(X^{'} Σ^{- 1} X)}^{- 1}, \end{matrix}

cov ({\hat{\hat{Θ}}}_{1}, {\hat{\hat{Θ}}}_{2}) = - {(X^{'} Σ^{- 1} X)}^{- 1} B_{1}^{'} Q_{1, 2},

var ({\hat{\hat{Θ}}}_{2}) = - Q_{2, 2},

and

(\begin{matrix} Q_{1, 1}, & Q_{1, 2} \\ Q_{2, 1}, & Q_{2, 2} \end{matrix}) = {(\begin{matrix} B_{1} {(X^{'} Σ^{- 1} X)}^{- 1} B_{1}^{'}, & B_{2} \\ B_{2}^{'}, & 0 \end{matrix})}^{- 1} .

Proof in [30].

3.4. Numerical Study

In the following example, the transformation parameters will be calculated in a simulated problem. We will consider a model of a building with four scans (coordinate systems) of measurement. See Figure 5 .

Let

Y_{1} = (\begin{matrix} Y_{1}^{I} \\ Y_{1}^{I I} \end{matrix})

Y_{2} = (\begin{matrix} Y_{2}^{I I} \\ Y_{2}^{I I I} \end{matrix})

Y_{3} = (\begin{matrix} Y_{3}^{I I I} \\ Y_{3}^{I V} \end{matrix})

Y_{4} = (\begin{matrix} Y_{4}^{I V} \\ Y_{4}^{I} \end{matrix})

For example,

Y_{1}^{I}

includes the measured coordinates on the I-side of an object that are given in the first stage of measurement. Vector

Y_{1}^{I I}

encompasses the measured coordinates on the II-side of an object that are given in the second stage of measurement.

Figure 5. Layout of measurement with 4 scans.

In Table 3, we present the measured coordinates:

We assume the coordinate accuracy of measured points in all scans to be approximately the same and to be given by a standard deviation of 10 mm.

For the above-mentioned algorithm, estimators of unknown parameters of translation and rotation will be calculated. During linearization, matrices obtained from partial derivations of the model according to individual parameters are used. The calculated matrices are given in Appendix A, Table A1.

Estimates computed by Formula (24) are

\begin{matrix} {\hat{T}}_{2} & = & (\begin{matrix} 0.55919 & 0.82904 & 0 \\ - 0.82904 & 0.55919 & 0 \\ 0 & 0 & 0 \end{matrix}), \\ {\hat{T}}_{3} & = & (\begin{matrix} - 0.91355 & 0.40674 & 0 \\ - 0.40674 & - 0.91355 & 0 \\ 0 & 0 & 0 \end{matrix}), \\ {\hat{T}}_{4} & = & (\begin{matrix} 0.46947 & - 0.88295 & 0 \\ 0.88295 & 0.46947 & 0 \\ 0 & 0 & 0 \end{matrix}), \\ {\hat{γ}}_{2} & = & {(99.3732, 10.8637, 0.0000)}^{'}, \\ {\hat{γ}}_{3} & = & {(89.6260, 84.0564, 9.9873)}^{'}, \\ {\hat{γ}}_{4} & = & {(2.5764, 95.9541, 15.0178)}^{'} . \end{matrix}

Furthermore, the resulting estimators of coordinates obtained in the system of scan 1 and absolute point positioning errors obtained from (25) are shown in Table 4.

4. Concluding Remarks

The paper attempts to present a state-of-the-art solution to the registration problem. Indeed, there are data registration problems that require high complexity for the laser systems and algorithms used. Many previously presented methods lack variability estimations. Our statistical model can serve as a complement to any registration algorithm since it provides an estimate of the variance matrix based on the knowledge of identical points. The method requires covariance matrix knowledge for identical point coordinates. The solution includes translational parameter estimation and the rotation matrix between the several coordinate systems. Unknown transformation parameters and the variance matrix were estimated by using the nonlinear regression model with nonlinear constraints. Separately, this approach is applicable in the case of geodetic measurements, when identical points are marked with reflective targets.

Author Contributions

Conceptualization, J.M. and P.C.; methodology, J.M. and P.C.; validation, , J.M. and P.C.; formal analysis, P.C.; investigation, P.C. and J.M.; writing-original draft preparation, P.C. and J.M.; supervision, P.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported by the Internal Grant Agency of the University of Pardubice. Vouchers from the review activity were also used to finance the contribution.

Informed Consent Statement

Not applicable.

Data Availability Statement

Some data are presented in tables in the article. Additional data sets are available upon request.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Table A1. Matrices

B_{1}

and

B_{2}

Table A1. Matrices

B_{1}

and

B_{2}

${\{B_{1}\}}_{i = 1}^{3 n_{2}} =$	$\frac{\partial}{\partial Θ_{1}^{'}} {[α_{2}^{(2)} - 1_{n_{2}} \otimes γ_{2} - (I_{n_{2}} \otimes T_{2}) β_{2}]}_{Θ_{1} = Θ_{1}^{(0)}, Θ_{2} = Θ_{2}^{(0)}} =$	$(0_{3 n_{2}, 3 n_{1}}, - (I_{n_{2}} \otimes T_{2}^{(0)}), I_{3 n_{2}}, 0_{3 n_{2}, 3 n_{3}}, 0_{3 n_{2}, 3 n_{3}}, 0_{3 n_{2}, 3 n_{4}}, 0_{3 n_{2}, 3 n_{4}}, 0_{3 n_{2}, 3 n_{1}}),$
${\{B_{1}\}}_{i = 3 n_{2} + 1}^{3 n_{2} + 3 n_{3}} =$	$\frac{\partial}{\partial Θ_{1}^{'}} {[α_{3}^{(2)} - 1_{n_{3}} \otimes γ_{2} - (I_{n_{3}} \otimes T_{2}) β_{3}]}_{Θ_{1} = Θ_{1}^{(0)}, Θ_{2} = Θ_{2}^{(0)}} =$	$(0_{3 n_{3}, 3 n_{1}}, 0_{3 n_{3}, 3 n_{2}}, 0_{3 n_{3}, 3 n_{2}}, I_{3 n_{3}}, 0_{3 n_{3}, 3 n_{3}}, 0_{3 n_{3}, 3 n_{4}}, 0_{3 n_{3}, 3 n_{4}}, 0_{3 n_{3}, 3 n_{1}}),$
${\{B_{1}\}}_{i = 3 n_{2} + 3 n_{3} + 1}^{3 n_{2} + 3 n_{3} + 3 n_{3}} =$	$\frac{\partial}{\partial Θ_{1}^{'}} {[α_{3}^{(3)} - 1_{n_{3}} \otimes γ_{3} - (I_{n_{3}} \otimes T_{3}) β_{3}]}_{Θ_{1} = Θ_{1}^{(0)}, Θ_{2} = Θ_{2}^{(0)}} =$	$(0_{3 n_{3}, 3 n_{1}}, 0_{3 n_{3}, 3 n_{2}}, 0_{3 n_{3}, 3 n_{2}}, 0_{3 n_{3}, 3 n_{3}}, I_{3 n_{3}}, 0_{3 n_{3}, 3 n_{4}}, 0_{3 n_{3}, 3 n_{4}}, 0_{3 n_{3}, 3 n_{1}}),$
${\{B_{1}\}}_{i = 3 n_{2} + 6 n_{3} + 1}^{3 n_{2} + 6 n_{3} + 3 n_{4}} =$	$\frac{\partial}{\partial Θ_{1}^{'}} {[α_{4}^{(3)} - 1_{n_{4}} \otimes γ_{3} - (I_{n_{4}} \otimes T_{3}) β_{4}]}_{Θ_{1} = Θ_{1}^{(0)}, Θ_{2} = Θ_{2}^{(0)}} =$	$(0_{3 n_{4}, 3 n_{1}}, 0_{3 n_{4}, 3 n_{2}}, 0_{3 n_{4}, 3 n_{2}}, 0_{3 n_{4}, 3 n_{3}}, 0_{3 n_{4}, 3 n_{3}}, I_{3 n_{4}}, 0_{3 n_{4}, 3 n_{4}}, 0_{3 n_{4}, 3 n_{1}}),$
${\{B_{1}\}}_{i = 3 n_{2} + 6 n_{3} + 3 n_{4} + 1}^{3 n_{2} + 6 n_{3} + 6 n_{4}} =$	$\frac{\partial}{\partial Θ_{1}^{'}} {[α_{4}^{(4)} - 1_{n_{4}} \otimes γ_{4} - (I_{n_{4}} \otimes T_{4}) β_{4}]}_{Θ_{1} = Θ_{1}^{(0)}, Θ_{2} = Θ_{2}^{(0)}} =$	$(0_{3 n_{4}, 3 n_{1}}, 0_{3 n_{4}, 3 n_{2}}, 0_{3 n_{4}, 3 n_{2}}, 0_{3 n_{4}, 3 n_{3}}, 0_{3 n_{4}, 3 n_{3}}, 0_{3 n_{4}, 3 n_{4}}, I_{3 n_{4}}, 0_{3 n_{4}, 3 n_{1}}),$
${\{B_{1}\}}_{i = 3 n_{2} + 6 n_{3} + 6 n_{4} + 1}^{3 n_{2} + 6 n_{3} + 6 n_{4} + 3 n_{1}} =$	$\frac{\partial}{\partial Θ_{1}^{'}} {[α_{1}^{(4)} - 1_{n_{1}} \otimes γ_{4} - (I_{n_{1}} \otimes T_{4}) β_{1}]}_{Θ_{1} = Θ_{1}^{(0)}, Θ_{2} = Θ_{2}^{(0)}} =$	$(- (I_{n_{1}} \otimes T_{4}), 0_{3 n_{1}, 3 n_{2}}, 0_{3 n_{1}, 3 n_{2}}, 0_{3 n_{1}, 3 n_{3}}, 0_{3 n_{1}, 3 n_{3}}, 0_{3 n_{1}, 3 n_{4}}, 0_{3 n_{1}, 3 n_{4}}, I_{3 n_{1}}),$
${\{B_{1}\}}_{i = 3 n_{1} + 3 n_{2} + 6 n_{3} + 6 n_{4} + 1}^{3 (n_{1} + n_{2}) + 6 (n_{3} + n_{4}) + 18} =$		$0_{18, 6 (n_{1} + n_{2} + n_{3} + n_{4})} .$
${\{B_{2}\}}_{i = 1}^{3 n_{2}} =$	$\frac{\partial}{\partial Θ_{2}^{'}} {[α_{2}^{(2)} - 1_{n_{2}} \otimes γ_{2} - (I_{n_{2}} \otimes T_{2}) β_{2}]}_{Θ_{1} = Θ_{1}^{(0)}, Θ_{2} = Θ_{2}^{(0)}} =$	$(0_{3 n_{2}, 3 n_{3}}, 0_{3 n_{2}, 3 n_{4}}, - 1_{n_{2}} \otimes I_{3}, 0_{3 n_{2}, 3}, 0_{3 n_{2}, 3}, - β_{2}^{'} \otimes I_{n_{3}}, 0_{3 n_{2}, 9}, 0_{3 n_{2}, 9}),$
${\{B_{2}\}}_{i = 3 n_{2} + 1}^{3 n_{2} + 3 n_{3}} =$	$\frac{\partial}{\partial Θ_{2}^{'}} {[α_{3}^{(2)} - 1_{n_{3}} \otimes γ_{2} - (I_{n_{3}} \otimes T_{2}) β_{3}]}_{Θ_{1} = Θ_{1}^{(0)}, Θ_{2} = Θ_{2}^{(0)}} =$	$((- I_{n_{3}} \otimes T_{2}), 0_{3 n_{3}, 3 n_{4}}, 0_{3 n_{3}, 3}, - 1_{n_{3}} \otimes I_{3}, 0_{3 n_{3}, 3}, - β_{3}^{'} \otimes I_{n_{3}}, 0_{3 n_{3}, 9}, 0_{3 n_{3}, 9}),$
${\{B_{2}\}}_{i = 3 n_{2} + 3 n_{3} + 1}^{3 n_{2} + 6 n_{3}} =$	$\frac{\partial}{\partial Θ_{2}^{'}} {[α_{3}^{(3)} - 1_{n_{3}} \otimes γ_{3} - (I_{n_{3}} \otimes T_{3}) β_{3}]}_{Θ_{1} = Θ_{1}^{(0)}, Θ_{2} = Θ_{2}^{(0)}} =$	$((- I_{n_{3}} \otimes T_{3}), 0_{3 n_{3}, 3 n_{4}}, 0_{3 n_{3}, 3}, - 1_{n_{3}} \otimes I_{3}, 0_{3 n_{3}, 3}, 0_{3 n_{3}, 9}, - β_{3}^{'} \otimes I_{n_{3}}, 0_{3 n_{3}, 9}),$
${\{B_{2}\}}_{i = 3 n_{2} + 6 n_{3} + 1}^{3 n_{2} + 6 n_{3} + 3 n_{4}} =$	$\frac{\partial}{\partial Θ_{2}^{'}} {[α_{4}^{(3)} - 1_{n_{4}} \otimes γ_{3} - (I_{n_{4}} \otimes T_{3}) β_{4}]}_{Θ_{1} = Θ_{1}^{(0)}, Θ_{2} = Θ_{2}^{(0)}} =$	$(0_{3 n_{4}, 3 n_{3}}, (- I_{n_{4}} \otimes T_{3}), 0_{3 n_{4}, 3}, - 1_{n_{4}} \otimes I_{3}, 0_{3 n_{4}, 3}, 0_{3 n_{4}, 9}, - β_{4}^{'} \otimes I_{n_{3}}, 0_{3 n_{4}, 9}),$
${\{B_{2}\}}_{i = 3 n_{2} + 6 n_{3} + 3 n_{4} + 1}^{3 n_{2} + 6 n_{3} + 6 n_{4}} =$	$\frac{\partial}{\partial Θ_{2}^{'}} {[α_{4}^{(4)} - 1_{n_{4}} \otimes γ_{4} - (I_{n_{4}} \otimes T_{4}) β_{4}]}_{Θ_{1} = Θ_{1}^{(0)}, Θ_{2} = Θ_{2}^{(0)}} =$	$(0_{3 n_{4}, 3 n_{3}}, (- I_{n_{4}} \otimes T_{4}), 0_{3 n_{4}, 3}, 0_{3 n_{4}, 3}, - 1_{n_{4}} \otimes I_{3}, 0_{3 n_{4}, 9}, 0_{3 n_{4}, 9}, - β_{4}^{'} \otimes I_{n_{3}}),$
${\{B_{2}\}}_{i = 3 n_{2} + 6 (n_{3} + n_{4}) + 1}^{3 n_{2} + 6 (n_{3} + n_{4}) + 3 n_{1}} =$	$\frac{\partial}{\partial Θ_{2}^{'}} {[α_{1}^{(4)} - 1_{n_{1}} \otimes γ_{4} - (I_{n_{1}} \otimes T_{4}) β_{1}]}_{Θ_{1} = Θ_{1}^{(0)}, Θ_{2} = Θ_{2}^{(0)}} =$	$(0_{3 n_{1}, 3 n_{3}}, 0_{3 n_{1}, 3 n_{4}}, 0_{3 n_{1}, 3}, 0_{3 n_{1}, 3}, - 1_{n_{1}} \otimes I_{3}, 0_{3 n_{1}, 9}, 0_{3 n_{1}, 9}, - β_{1}^{'} \otimes I_{3})$ .

References

Goshtasby, A.A. Image Registration; Advances in Pattern Recognition; Springer: London, UK, 2012. [Google Scholar] [CrossRef]
Pomerleau, F.; Colas, F.; Siegwart, R. A Review of Point Cloud Registration Algorithms for Mobile Robotics. Found. Trends Robot. 2015, 4, 1–104. [Google Scholar] [CrossRef]
Cheng, L.; Chen, S.; Liu, X.; Xu, H.; Wu, Y.; Li, M.; Chen, Y. Registration of laser scanning point clouds: A review. Sensors. Multidisciplinary Digital Publishing Institute. Sensors 2018, 18, 1641. [Google Scholar] [CrossRef]
Dong, Z.; Liang, F.; Yang, B.; Xu, Y.; Zang, Y.; Li, J.; Wang, Y.; Dai, W.; Fan, H.; Hyypp, J.; et al. Registration of large-scale terrestrial laser scanner point clouds: A review and benchmark. ISPRS J. Photogramm. Remote Sens. 2020, 163, 327–342. [Google Scholar] [CrossRef]
Kopacik, A.; Erdeeyi, J.; Kyrinovic, P. Engineering Surveys for Industry, Chapter Terrestrial Laser Scanning Systems; Springer: Berlin/Heidelberg, Germany, 2020; ISBN 9783030483081. [Google Scholar] [CrossRef]
Chmelar, P.; Rejfek, L.; Nguyen, T.N.; Ha, D. Advanced Methods for Point Cloud Processing and Simplification. Appl. Sci. 2020, 10, 3340. [Google Scholar] [CrossRef]
Chen, Y.; Gerard, M. Object modelling by registration of multiple range images. Image Vision Comput. 1991, 10, 145–155. [Google Scholar] [CrossRef]
Besl, P.J.; McKay, N.D. A Method for Registration of 3-D Shapes. IEEE Trans. Pattern Anal. Mach. Intell. 1992, 14, 239–256. [Google Scholar] [CrossRef]
Nüchter, A. 3D Robotic Mapping: The Simultaneous Localization and Mapping Problem with Six Degrees of Freedom; Springer: Berlin/Heidelberg, Germany, 2009; ISBN 978-354089883-2. ISSN 16107438. [Google Scholar] [CrossRef]
Nüchter, A.; Elseberg, J.; Schneider, P.; Paulus, D. Study of parameterizations for the rigid body transformations of the scan registration problem. Comput. Vis. Image Underst. 2010, 114, 963–980. [Google Scholar] [CrossRef]
Peter, B.; Wolfgang, S. The normal distributions transform: A new approach to laser scan matching. In Proceedings of the 2003 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS 2003) (Cat. No. 03CH37453), Las Vegas, NV, USA, 27–31 October 2003; Volume 3, pp. 2743–2748. [Google Scholar]
Magnusson, M.; Lilienthal, A.; Duckett, T. Scan registration for autonomous mining vehicles using 3D-NDT. J. Field Robot. 2007, 24, 803–827. [Google Scholar] [CrossRef]
Magnusson, M. The Three-Dimensional Normal-Distributions Transform an Efficient Representation for Registration, Surface Analysis, and Loop Detection. In Studies in Technology; Örebro University: Örebro, Sweden, 2013. [Google Scholar]
Beran, L.; Chmelar, P.; Rejfek, L. Navigation of Robotics Platform Using Advanced Image Processing Navigation Methods. In Proceedings of the VIPIMAGE (V. Eccomas Thematic Conference on Computional Vision and Medical Image Processing), Tenerife, Spain, 19–21 October 2015; ISBN 978-113802926-2. [Google Scholar]
Beran, L.; Chmelar, P.; Rejfek, L. Navigation of Robotics Platform using Monocular Visual Odometry. In Proceedings of the Radioelektronika 25th International Conference, Pardubice, Czech Republic, 21–22 April 2015. [Google Scholar]
Krajnik, T. Large-Scale Mobile Robot Navigation and Map Building. Ph.D. Thesis, Czech Technical University in Prague, Faculty of Electrical Engineering, Department of Cybernetics, Prague, Czech Republic, 2011. [Google Scholar]
Lu, F.; Milios, E. Robot pose estimation in unknown environments by matching 2D range scans. J. Intell. Robot. Syst. 1997, 18, 249–275. [Google Scholar] [CrossRef]
Montesano, L.; Minguez, J.; Montano, L. Probabilistic scan matching for motion estimation in unstructured environments. In Proceedings of the IEEE International Conference on Intelligent Robots and Systems (IROS), Kyoto, Japan, 23–27 October 2005. [Google Scholar] [CrossRef]
Hähnel, D.; Burgard, W. Probabilistic Matching for 3D Scan Registration. In Proceedings of the VDI-Conference Robotik 2002 (Robotik), Forum am Schlosspark, Germany, 19–20 June 2002; Available online: https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=4f75beeeb9958a948492490f69c25f37165a3908 (accessed on 10 August 2023).
Boughobel, F.; Muharrem, M.; Koschan, A.; Abidi, M. A new method for the registration of three-dimensional point-sets: The Gaussian Fields framework. Image Vis. Comput. 2010, 28, 124–137. [Google Scholar] [CrossRef]
Boughorbel, F.; Koschan, A.; Abidi, B.; Abidi, M. Gaussian fields: A new criterion for 3D rigid registration. Pattern Recognit. 2004, 37, 1567–1571. [Google Scholar] [CrossRef]
Mitra, N.J.; Gelfand, N.; Pottmann, H.; Guibas, L. Registration of point cloud data from a geometric optimization perspective. In Proceedings of the Symposium on Geometry Processing, Nice, France, 8–10 July 2004; pp. 22–31. [Google Scholar] [CrossRef]
Burguera, A.; Gonz lez, Y.; Oliver, G. The likelihood field approach to sonar scan matching. In Proceedings of the IEEE International Conference on Intelligent Robots and Systems (IROS), Kyoto, Japan, 23–27 October 2008; pp. 2977–2982. [Google Scholar] [CrossRef]
Lafferty, J.; McCallum, A.; Pereira, F. Conditional randomfields: Probabilistic models for segmenting and labeling sequence data. In Proceedings of the Eighteenth International Conference on Machine Learning (ICML), San Francisco, CA, USA, 28 June–1 July 2001; ISBN 1-55860-778-1. [Google Scholar]
Forsman, P.; Halme, A. Feature based registration of range images for mapping of natural outdoor environments. In Proceedings of the International Symposium on 3D Data Processing, Visualization, and Transmission (3DPVT), Washington, DC, USA, 14–16 June 2004. [Google Scholar] [CrossRef]
Johnson, A.E. Spin Images: A Representation for 3-D Surface Matching. Ph.D. Thesis, Carnegie Mellon University, Pittsburgh, PA, USA, 1997. [Google Scholar]
Olsen, M.; Johnstone, E.; Kuester, F.; Driscoll, N.; Ashford, S. New Automated Point-Cloud Alignment for Ground-Based Light Detection and Ranging Data of Long Coastal Sections. J. Surv. Eng. 2011, 137, 14–25. [Google Scholar] [CrossRef]
Ochmann, S.; Vock, R.; Wessel, R.; Klein, R. Automatic reconstruction of parametric building models from indoor point clouds. Comput. Graph. 2016, 54, 94–103. [Google Scholar] [CrossRef]
Marek, J.; Rak, J.; Jetensky, P. Statistical solution of 3D transformation problem. In Computer Science Research Notes, CSRN 2503; Zapadoceska Univerzita v Plzni: Plzen, Czech Republic, 2015; pp. 85–90. ISBN 978-808694367-1. [Google Scholar]
Kubacek, L. Statistical Theory of Geodetic Networks; VUGTK: Zdiby, Czech Republic, 2013. [Google Scholar]

Figure 1. Measuring device.

Figure 2. One measurement frames of Figure 1.

Figure 3. 3D point cloud from 1st and 2nd position.

Figure 4. ICP algorithm result. Left side the registered point clouds from Figure 3 by different colors. Right side the colored version.

Table 1. Basic methods comparison part 1.

Algorithm	Based on (Points/Model)	Covarience Matrix	Important Requirements	Speed Processing	Computation Complexity
3D ICP	points	is estimated	-	slow for large datasets	$O (N^{2})$
3D NDT	3D normal distribution	is estimated	-	faster than ICP	$O (N)$
Feature based	search space features	depends on application	varied environment	fast—depends on selected feature	depends on selected descriptor
IDC	points	is estimated	-	faster than ICP	similar to ICP but not specified
pIC	points	is estimated	-	slow for large datasets	similar to ICP but not specified
Point-based probabilistic	Probability functions	ss estimated	-	faster than ICP	similar to ICP but not specified
Gaussian fields	Gaussian mixture	is estimated	-	faster than ICP	$O (N)$
Quadratic patches	Quadratic approximation	is estimated	-	faster than ICP	not specified
Likelihood-fields	Likelihood-fields	uses identity matrix	-	similar to NDT	not specified
CRF matching	Conditional random fields	uses log-likelihood function	make sense only in 2D	slow for large datasets	$O (N)$ in current scan, $O (N^{2})$ in reference scan
Branch and bound	Rotation symmetric features	is necessary, get from sensors	make sense only in 2D	not specified	not specified
Local geometric features	Local geometric features	Is estimated	-	Depending on used algorithm	Depending on used algorithm
PointReg	points	is necessary	only consistent object, moving has to be eliminated	slow due to user interaction	not specified

Table 2. Basic methods comparison part 2.

Algorithm	Precision	Advantages	Necessity of Human Intervention	Suitable Application	Origin
3D ICP	good - depens on input data	precise with suitable data	initial parameters setup	precise 3D point clouds	1992
3D NDT	good with suitable settings	space modeling, using of less accurate data	initial parameters setup	various input data	2007
Feature based	good, depens on environment	working with camera images	initial parameters setup	robotic navigation	1998, 2006 (real time)
IDC	more robust, less precision than ICP	robustnes for bigger initial pose	initial parameters setup	similar to ICP	1994
pIC	better accuracy, robustness and convergence than ICP and IDC	incorporates sensors quality and noise	initial parameters setup	similar to ICP	2005
Point-based probabil	better than ICP	scans treating as probability function	initial parameters setup	large 3D scans	2002
Gaussian fields	good with suitable settings	similar to 3D NDT	initial parameters setup	various input data	2004
Quadratic patches	good with suitable settings	similar to 3D NDT	initial parameters setup	various input data	2004
Likelihood-fields	similar to NDT	possiblity to extend scans to get more points	initial parameters setup	various input data	2008
CRF matching	good with suitable settings	using od different user-features, robust to initial pose error	initial parameters setup	only for 2D scans	2001
Branch and bound	not specifed, dependent on set of common symetric features	using on highly unstructured environments	initial parameters setup	natural outdoor environment	2004
Local geometric feat	good depending on input data	matching local surfaces	initial parameters setup	point clouds wit varied objects	1997
PointReg	low RMS with the suitable settings and data handling	good with good data handling	manual	pitfalls 3D range scans	2011

Table 3. Measurements.

Y $_{1}^{I}$	Y $_{1}^{II}$	Y $_{2}^{II}$	Y $_{2}^{III}$	Y $_{3}^{III}$	Y $_{3}^{IV}$	Y $_{4}^{IV}$	Y $_{4}^{I}$
629.1793	629.0871	521.6030	533.7651	−457.5236	−446.6342	208.9794	223.0744
84.1583	84.0498	−464.5526	−469.4856	−258.9134	−265.6701	702.7168	691.0423
15.1976	4.0980	4.0980	16.6208	26.6208	17.8997	22.8997	30.1976
626.6654	638.7188	532.2563	536.3855	−451.5327	−445.1702	209.2158	218.5299
87.9686	90.5740	−468.8736	−462.9402	−262.6306	−263.4124	700.0364	690.6115
3.5058	7.8426	7.8426	7.7437	17.7437	12.0001	17.0001	18.5058
622.4489	640.0699	533.7651	536.7300	−450.7452	−444.5112	209.3221	210.9074
94.3597	91.4826	−469.4856	−462.0798	−263.1193	−262.3962	698.8299	689.8890
4.0666	5.5998	5.5998	3.3267	13.3267	14.6452	19.6452	19.0666
622.0117			536.5612	−451.1310	−443.3154	209.5151	210.1170
95.0224			−462.5013	−262.8799	−260.5522	696.6406	689.8140
3.0727			16.2879	26.2879	17.3570	22.3570	18.0727
			538.5282	−446.6342	−445.0099	209.2416
			−457.5883	−265.6701	−263.1652	699.7429
			17.6355	27.6355	20.7081	25.7081
					−439.5869	210.1170
					−254.8022	689.8140
					14.7198	19.7198
$n_{1} = 4$	$n_{2} = 3$		$n_{3} = 5$		$n_{4} = 6$		$n_{1} = 4$

Table 4. Estimates of coordinates.

${\hat{β}}_{1}$	$^{s}$ ${\hat{β}}_{1}$	${\hat{β}}_{2}$	$^{s}$ ${\hat{β}}_{2}$	${\hat{β}}_{3}$	$^{s}$ ${\hat{β}}_{3}$	${\hat{β}}_{4}$	$^{s}$ ${\hat{β}}_{4}$
629.2236	0.0134	629.2727	0.0151	639.9965	0.0254	632.7031	0.1205
84.0842	0.0162	84.1765	0.0108	91.4590	0.0235	102.1926	0.1057
14.8946	0.0180	4.2584	0.0149	16.7208	0.0232	7.9921	0.1125
626.7046	0.0121	638.7539	0.0133	636.1058	0.0753	630.6615	0.0852
87.9179	0.0107	90.4175	0.0154	97.5176	0.0541	100.7179	0.0816
3.4601	0.0155	7.8524	0.0151	7.5773	0.0429	2.0002	0.0344
622.3238	0.0142	640.0670	0.0180	635.6114	0.0907	629.4778	0.0233
94.3276	0.0139	91.4741	0.0189	97.8561	0.0687	99.8918	0.0422
4.1908	0.0188	5.6040	0.0201	3.2677	0.0789	4.6397	0.0179
621.9169	0.0190			635.8859	0.0595	627.7632	0.0912
95.0236	0.0224			97.9292	0.0614	98.7844	0.0894
2.9661	0.0256			16.2601	0.0209	7.4481	0.0796
				632.8328	0.1003	630.4326	0.0488
				102.1181	0.0982	100.3622	0.0198
				17.6777	0.1143	10.7676	0.0158
						622.0444	0.0129
						94.8902	0.0164
						4.7548	0.0246

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2023 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 (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Marek, J.; Chmelař, P. Survey of Point Cloud Registration Methods and New Statistical Approach. Mathematics 2023, 11, 3564. https://doi.org/10.3390/math11163564

AMA Style

Marek J, Chmelař P. Survey of Point Cloud Registration Methods and New Statistical Approach. Mathematics. 2023; 11(16):3564. https://doi.org/10.3390/math11163564

Chicago/Turabian Style

Marek, Jaroslav, and Pavel Chmelař. 2023. "Survey of Point Cloud Registration Methods and New Statistical Approach" Mathematics 11, no. 16: 3564. https://doi.org/10.3390/math11163564

APA Style

Marek, J., & Chmelař, P. (2023). Survey of Point Cloud Registration Methods and New Statistical Approach. Mathematics, 11(16), 3564. https://doi.org/10.3390/math11163564

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Survey of Point Cloud Registration Methods and New Statistical Approach

Abstract

1. Introduction

1.1. Optical Rangefinder

1.2. 3D Point Cloud

2. A Survey of Methods

2.1. ICP Algorithm

2.1.1. ICP Algorithm Results

2.1.2. ICP Algorithm Problem

2.2. Normal Distribution Transform (NDT)

2.3. Feature Based Registration

2.4. IDC

2.5. pIC

2.6. Point-Based Probabilistic Registration

2.7. Gaussian Fields

2.8. Quadratic Patches

2.9. Likelihood-Field Matching

2.10. CRF Matching

2.11. Branch-and-Bound Registration

2.12. Registration Using Local Geometric Features

2.13. PoitnReg

2.14. The Automatic Mapping of Parametric Objects (Buildings) from Indoor Point Clouds

2.15. Basic Registration Methods Comparison

3. The Optimal Statistical Transformation Model Taking into Account Errors in Coordinate Systems

3.1. Model with Four Scans

3.2. Linearization of a Model with 4 Scans

3.3. The Estimates of Unknown Parameters in a Model with 4 Scans

3.4. Numerical Study

4. Concluding Remarks

Author Contributions

Funding

Informed Consent Statement

Data Availability Statement

Conflicts of Interest

Appendix A

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI