Composite Ski-Resort Registration Method Based on Laser Point Cloud Information

Abstract

The environment of ski resorts is usually complex and changeable, and there are few characteristic objects in the background, which creates many difficulties for the registration of ski-resort point cloud datasets. However, in the traditional iterative closest point (ICP) algorithm, two points need to have good initial positions, otherwise it is easy to get caught up in local optimizations in registration. Aiming at this problem, according to the topographic features of ski resorts, this paper put forward a ski-resort coarse registration method based on extraction, and matching between feature points is proposed to adjust the initial position of two point clouds. Firstly, the feature points of the common part of the point cloud datasets are extracted based on the SIFT algorithm; secondly, the Euclidean distance between the feature normal vectors is used as the pairing condition to complete the pairing between the feature points in the point cloud datasets; then, the feature point pair is purified by using the included angle of the normal vector; finally, in the process of coarse registration, the rotation matrix and translation vector between point clouds are solved by the unit quaternion method. Experiments demonstrate that the proposed coarse registration method based on the normal vector of feature points is helpful to the smooth completion of the subsequent fine registration process, avoids the phenomenon of falling into local optimization, and effectively completes the ski-resort point cloud registration.

1. Introduction

With the successful holding of the 24th Beijing Winter Olympic Games in 2022, the enthusiasm of the masses to participate in winter sports has been greatly stimulated and, for all we know, the more typical winter sport is skiing. Accordingly, scientific research and participants related to skiing have achieved rapid growth in varying degrees [1,2,3].

In recent years, because of its high accuracy and efficiency, three-dimensional laser radar (Light Detection and Ranging, LiDAR) has been widely used in computer vision, cultural relic protection, reverse engineering and geographic mapping [4,5,6,7,8]. Additionally, according to the many advantages of LiDAR, it has become a new application direction to use the UAV LiDAR system to collect high-precision three-dimensional datasets of ski resorts [9,10], complete three-dimensional real scene reconstruction and realize snow track information extraction. The acquisition of 3D digital ski resorts is conducive to the daily maintenance of ski resorts (such as snow making, snow pressing, route planning, etc.) and the generation of e thvirtual scene in simulated training and entertainment system [11], which has great application prospects. The effective registration of ski-resort point cloud datasets collected by different sorties and routes has become an important foundation and key step of the ski-resort 3D reconstruction. Besides, because of the narrow visual surface of LiDAR sensor and the high spatial complexity of geographical entities, data acquisition based on LiDAR technology usually needs to be carried out along multiple perspectives [12,13], and then accomplishes the comprehensive expression of geographical entities based on the corresponding registration algorithm.

According to the different feature search space, the registration process can be summarized as global registration and local registration [14,15,16,17]. Global registration is usually based on local geometric features that can be uniquely matched without assuming an initial position. These features can be prearranged targets, scanning corners, lines, or planes in the scene [18]. The registration efficiency is high, but if the geometric conditions are insufficient and the feature distribution is uneven, the reliability of the registration results will be reduced. Local registration is to use the original point cloud to optimize the transformation matrix to obtain high matching accuracy when the position estimation of the two-point set is known [19]. Among them, the ICP algorithm is one of the most commonly used registration algorithms because of its simple steps, high precision and good robustness [20,21]. It has high requirements for the initial position between point clouds, and when the initial value error is large, it is easy to produce local optimization. Because the ICP algorithm directly uses the nearest sampling point in the point cloud datasets to establish the corresponding relationship, it has high requirements for the overlap of point clouds and takes a long time to search, so it cannot meet the application requirements when the amount of data is large [22,23,24,25].

Over the years, a considerable number of researchers and engineers at home and abroad have made a large number of improvements to the ICP algorithm. Men et al. changed the selection method of corresponding points by the Euclidean distance of point to point to the Euclidean distance of point to tangent plane, which improved the accuracy of registration, but the efficiency is still not high [26]. Koide et al. proposed the iterative dual correspondence (IDC) algorithm, which integrates the correspondence rules of the nearest point and matching distance point [27]. When the initial position is far away, IDC has better convergence than ICP, but when the initial position is close, its accuracy is lower than the ICP algorithm. Lei et al. put forward the algorithm of iterating the nearest segment and the nearest patch on the basis of ICP [28]. Firstly, the points in the two datasets are directly connected or triangulated, and then the corresponding line segment or triangular patch in the two datasets is found according to certain criteria, and the objective equation is established to calculate the rotation matrix. However, the stability of the algorithm used to find the corresponding relation criterion is yet to be verified. Rusu controlled the search range of the corresponding point set in the pre-divided cube grid, not only considering the curvature change of the point cloud, but also taking the information such as intensity, color and texture as the standard of point cloud registration [29].

In conclusion, considering that most registration algorithms converge on the assumption that the initial values of registration parameters are known and good, and when the initial values are incorrect, the algorithm cannot get the correct results; therefore, how to determine the initial value of the registration is a key. This paper presents a ski-resort point cloud registration algorithm based on the Euclidean distance and the included angle between the normal vectors of feature points. In the coarse registration stage, firstly, SIFT algorithm is used to extract feature points. Secondly, the Euclidean distance between the normal vectors of feature points is used as the judgment basis to automatically match the feature points of the ski-resort point cloud datasets. Then, the feature points are purified according to the topological relationship, that is, the included angle between the normal vectors of the feature points. Finally, the unit quaternion method is used to rotate and translate the two point clouds, so as to provide a good results for the subsequent fine registration work, which can effectively solve the problem of falling into local optimization and shorten the time of fine registration.

2. Study Area and Data

The LiDAR datasets used in this paper are acquired through the unmanned aircraft LiDAR System (RIEGL VUX-1 UAV System). The collection site is located in the Zhangjiakou Wanlong ski resort, the most famous ski resort in China, near Yanqing District, Beijing (249 km away from Beijing, 2110.3 m above sea level, more than 200,000 square meters). The collection time was 12:00 to 14:00 noon, and the temperature of the ski resort was −18 °C. Figure 1a shows the environmental conditions of the Wanlong ski resort at the time of collection, Figure 1b is an aerial picture taken from the air during acquisition, and Figure 2 is an introduction to the system architecture.

The seven ski-tracks of the Wanlong ski resort are all made of artificial snow, with high hardness. The environment of the snow field is relatively complex, including not only five main ski tracks, but also two training tracks, as well as trees, riprap, etc. DJI M600 PRO UAV is used as the flight platform with LiDAR, because of its excellent stability and flight performance, and its endurance time can reach 30 min. In the acquisition process, the scanning frequency of LiDAR was set to 550 KHz, the overlap rate of sidebands was 20 percent, the flight altitude was 110 m, the flight speed was 5 m/s and a total of two sorties were flown to complete the data acquisition.

3. Methodology

3.1. Coarse Registration of Point Cloud in Ski Resort

In order to get a good initial position, two pieces of the point cloud of the Wanlong ski resort need to be coarse-registered before fine registration.

3.1.1. SIFT Algorithm to Extract Feature Points

SIFT is a local feature descriptor proposed by David Lowe in 1999, and in 2004, has carried on further development and improvement [30,31,32,33]. SIFT features keep good invariance to rotation, translation, scale scaling, occlusion and brightness change, and are stable to visual change and affine transformation. In addition, it is especially suitable for high-speed and efficient matching in massive databases. The principle of SIFT algorithm can be simply described as: a picture in different scale space $L(x,y,\sigma )$ is defined as the convolution of a gauss function $G(x,y,\sigma )$ and the original picture $I(x,y)$. The relationship is as follows:

$$L(x,y,\sigma )=G(x,y,\sigma )\ast I(x,y)$$

(1)

In Formula (1), $(x,y)$ represents the pixel position of the picture; ∗ is a convolution operation; $\sigma$ is the scale space factor.An efficient gauss difference operator $D(x,y,\sigma )$ is used for extremum detection to select stable key points, as follows:

$$D(x,y,\sigma )=\left(G\right(x,y,k\sigma )-G(x,y,\sigma \left)\right)\ast I(x,y)=L(x,y,k\sigma )-L(x,y,\sigma )$$

(2)

In Formula (2), k is a constant. The adjacent layers in DoG space are detected, and some extremum points that are detected together constitute the feature points. By comparing each pixel with 8 points of the same scale and 18 points of the adjacent scale, the extreme points of the DoG function can be found. If the DoG operator value of the pixel is an extreme value in the 26 neighborhood points, this pixel is defined as a feature point. At present, the feature points are the extreme points in the discrete space. Through the fitting of three-dimensional surface, the position and scale of feature points are accurately determined, low contrast feature points and unstable edge points are removed, the matching stability is enhanced and the anti noise ability is improved.

The least square fitting is carried out by using the DoG function in Taylor quadratic expansion:

$$D\left(X\right)=D+\frac{\partial D^T}{\partial X}X+\frac{1}{2}{X}^T\frac{{\partial }^{2}D}{\partial X^2}$$

(3)

where $X={(x,y,\sigma )}^T$. If Equation (3) are equal to zero, we can get the offset of the extreme point:

$$\widehat{X}=-\frac{{\partial }^2{D}^{-1}}{\partial X^2}\frac{\partial D}{\partial X}$$

(4)

where $\widehat{X}={(x,y,\sigma )}^T$ is the relative offset of the interpolation center. If the offset of any dimension is not less than 0.5 (that is x, y, or $\sigma$), this indicates that the interpolation center is already at other adjacent points, and the current critical point must be changed. Therefore, $\widehat{X}$ is used to calibrate the coordinates of the original feature points, and then the new coordinates after position correction are obtained.

3.1.2. Feature Point Matching of Point Cloud Datasets

After the feature points of the common part of the two point clouds are extracted, the feature points would be matched, and the Euclidean distance between the normal vectors is used to judge the similarity of the feature points in the cloud. Then, the Euclidean distance between the normal vector $n_i$ of the point cloud and $n_j$ of the feature points of the target cloud is calculated, and the two points with the smallest Euclidean distance are selected by comparing each other. If the preset threshold $\tau (\tau >0)$ is greater than the minimum distance, then the nearest point pairs of the normal Euclidean distance between the feature points in the target and source point clouds are feature pairs. If the threshold is reduced, the number of matching point pairs will decrease, but the stability of matching will increase.

It should be noted that the normal vector of each feature points can be divided into positive and negative. Therefore, In order to keep the direction of the normal vector of the scanned object surface consistent, all the normal vectors are adjusted to satisfy $n_i·n_j<0(i\ne j)$, so that the matching of feature points is more accurate.

The solution of normal vectors are based on the surface of points, ski-resort point cloud datasets are discrete points, and the information recorded in point clouds is the 3D coordinates of each independent point. Therefore, in order to obtain the normal vector of each feature point, the points within a certain radius near the point are fitted to the surface, and the normal vector of the target point is obtained on the basis of the fitting surface. The calculation process is as follows:

(1) For the n points $(x_i,y_i,z_i),i=1,2,\cdots ,n$ to be fitted, the fitted plane equation is:

$$ax+by+cz=d,d⩾0$$

(5)

where $a,b,c$ should satisfy $a^2+b^2+c^2=1$. In addition, the distance from any point $(x_i,y_i,z_i)$ to the plane is:

$$d_i=|a{x}_i+b{y}_i+c{z}_i-d|$$

(6)

(2) In order to obtain the best fitting plane, to meet the needs of the least squares, it is necessary to obtain the minimum sum of squares of all distances between each point and the plane, ${\displaystyle e=\sum ^n{d}_i^2\to min}$, so that the problem of finding the normal vector can be transformed into that of finding the extremum:

$${\displaystyle f=\sum ^n{d}_i^2-\mu (a^2+b^2+c^2-1)}$$

(7)

(3) Taking the partial derivative of f with respect to four unknown parameters $d,a,b,c$, and then getting the unknown parameter d:

$$d=a\frac{{\displaystyle \sum _{i=1}^n{x}_i}}{n}+b\frac{{\displaystyle \sum _{i=1}^n{y}_i}}{n}+c\frac{{\displaystyle \sum _{i=1}^n{z}_i}}{n}$$

(8)

(4) The distance from any point to the plane can be rewritten as:

$$d=\left|a\left(x_i-\frac{{\displaystyle \sum _{i=1}^n{x}_i}}{n}\right)+b\left(y_i-\frac{{\displaystyle \sum _{i=1}^n{y}_i}}{n}\right)+c\left(z_i-\frac{{\displaystyle \sum _{i=1}^n{z}_i}}{n}\right)\right|$$

(9)

Let $\overline{x}=\frac{{\displaystyle \sum _{i=1}^n{x}_i}}{n},\overline{y}=\frac{{\displaystyle \sum _{i=1}^n{y}_i}}{n},\overline{z}=\frac{{\displaystyle \sum _{i=1}^n{z}_i}}{n},$ thereby:

$$d_i=\left|a(x_i-\overline{x})+b(y_i-\overline{y})+c(z_i-\overline{z})\right|$$

(10)

(5) Finding the partial derivative of formula (10) with respect to $a,b$ and c:

$${\displaystyle \frac{\partial f}{\partial a}=2\sum _{i=1}^n\left(a(x_i-\overline{x})+b(y_i-\overline{y})+c(z_i-\overline{z})\right)(x_i-\overline{x})-2\mu a}$$

(11)

Let $\Delta x_i=x_i-\overline{x},\Delta y_i=y_i-\overline{y},\Delta z_i=z_i-\overline{z}$, it can be transformed into:

$${\displaystyle \frac{\partial f}{\partial a}=2\sum _{i=1}^n\left(a\Delta x_i+b\Delta y_i+c\Delta z_i\right)\Delta x_i-2\mu a}$$

(12)

In addition,

$${\displaystyle \frac{\partial f}{\partial b}=2\sum _{i=1}^n\left(a\Delta x_i+b\Delta y_i+c\Delta z_i\right)\Delta x_i-2\mu b}$$

(13)

$${\displaystyle \frac{\partial f}{\partial c}=2\sum _{i=1}^n\left(a\Delta x_i+b\Delta y_i+c\Delta z_i\right)\Delta x_i-2\mu c}$$

(14)

(6) Combining Equation (12) to (14), it can be rewritten as:

$$\left[\begin{array}{ccc}\sum \Delta x_i\Delta x_i& \sum \Delta x_i\Delta y_i& \sum \Delta x_i\Delta z_i\\ \sum \Delta x_i\Delta y_i& \sum \Delta y_i\Delta y_i& \sum \Delta y_i\Delta z_i\\ \sum \Delta x_i\Delta z_i& \sum \Delta y_i\Delta z_i& \sum \Delta z_i\Delta z_i\end{array}\right]\left[\begin{array}{c}a\\ b\\ c\end{array}\right]=\mu \left[\begin{array}{c}a\\ b\\ c\end{array}\right]$$

(15)

Therefore, ${(a,b,c)}_T$ is the normal vector of the point in the plane.

If a feature point is selected from the target ski-resort datasets and its normal vector is ${(a_j,b_j,z_j)}^T$, the formula is the Euclidean distance between two normal vectors:

$$d_{ij}=\sqrt{{(a_i-a_j)}^2+{(b_i-b_j)}^2+{(c_i-c_j)}^2}$$

(16)

3.1.3. Purification of Feature Points

The threshold set in the feature point pairing cannot completely and correctly correspond to the feature points one by one, and the phenomenon of false matching cannot be avoided. In theory, as long as 3 pairs of correct matching point pairs are extracted, an accurate rotation and translation matrix of ski-resort point cloud datasets can be obtained. However, if one pair of feature points is mismatched, there will be a great difference between the calculated transformation matrix and the actual transformation matrix. Therefore, in order to further improve the accuracy of the rotation matrix and translation vector, it is necessary to further purify the paired feature point pairs.

Currently, the most common algorithm for matching point refinement is the RANSAC algorithm (random sample consensus, RANSAC), which iterates through a set of observations containing outliers to estimate the parameters of the mathematical model. In addition, the point cloud datasets can be divided into valid data and invalid data. The RANSAC algorithm builds different models based on different data, the data satisfies the model, it is valid data; otherwise, it is invalid data. The results demonstrate that the model data obtained by this algorithm is robust and the parameters estimated by this algorithm are relatively accurate. However, the model fitted by this algorithm is different for different data every time, so it is not suitable for complex point cloud scene. Moreover, no upper limit is set on the number of iterations of the computational model. If the number of iterations is limited, the result may not be optimal or even wrong. Therefore, this paper adopts a matching point purification algorithm based on the included angle of the normal vector between the feature point pairs as the evaluation quantity. The included angle of normal vector between adjacent point clouds is affected by the slope value of the surface. The greater the slope, the greater the included angle. On the contrary, the smaller the slope, the smaller the included angle. The specific relationship is shown in Figure 3. Although the source ski-resort datasets and the target ski-resort datasets are in different spatial coordinate systems, their spatial topological relationship should be consistent. In addition to translation changes, there should also be a rotation relationship between the corresponding point clouds, and the included angle between the normal vectors should meet a certain mathematical relationship. Therefore, a threshold G is set in advance for the angle between the feature points and the normal vector to eliminate the mismatched points, so as to achieve the effect of purification.

For any pair of feature points A and B, their normal vectors are $n_A$ and $n_B$, respectively. The greater the difference between the normal vectors, the smaller the cosine value of the included angle, and the smaller $n_A·n_B$. Therefore, the threshold G is set according to the cosine of the included angle between the normal vectors, and the point pairs whose product of the normal vectors is less than G are regarded as mismatched points and eliminated.

$$n_A·n_B⩾G$$

(17)

3.1.4. Point Cloud Coordinate Transformation

When the feature points of the source ski-resort datasets and the target ski-resort datasets are matched one by one, the translation and rotation matrix should be calculated by using the feature point pairs. The commonly used three-dimensional coordinate transformation method obtains seven parameters by using the Boolean model, including three coordinate translation parameters, three angular rotation parameters and one scale parameter. Because the point cloud registration does not involve the size transformation of the point cloud scale, only the first six parameters are used to translate and rotate the source point cloud. However, this method contains the unstable error after linearization, the error of conversion parameter can be ignored only in the case of small angle conversion, and the precision cannot meet the requirement of conversion when the coordinate conversion angle is large. Therefore, the unit quaternion method is selected in this paper, which is more suitable for large angle three-dimensional coordinate transformation.

Quaternions are proposed by 1843, whose essence is a simple hypercomplex consisting of 1 real part unit and 3 imaginary part units [34,35,36]. When using the unit quaternion method to convert the three-dimensional coordinates of the point cloud, the common part of the two point clouds must be found, the rotation matrix and the translation vector of the overlapping part of ski-resort point cloud datasets must be obtained, so as to be applied to the conversion of the whole part of the source ski-resort datasets. The specific conditions for its use are: (1) The common parts A and B of the ski-resort point cloud datasets are selected, respectively, and the number of point clouds in A and B are equal; (2) the point cloud of A and B should meet the one-to-one correspondence, which means that $A_i$ and $B_i$ represent the same point in different coordinate systems. In the former part, the feature points of the common part of the source ski-resort datasets and the target ski-resort datasets are matched one by one, and the matching point pair can meet the requirements of the unit quaternion method. The specific algorithm of the unit quaternion method is as follows:

(1) Calculation the center of gravity of A and B point clouds, respectively, and the center of gravity coordinates are $Z_A$ and $Z_B$.

$$\left\{\begin{array}{cc}\hfill {\displaystyle Z_A=\frac{1}{N}\sum _{i=1}^{N_A}{A}_i}& \\ \hfill {\displaystyle Z_B=\frac{1}{N}\sum _{i=1}^{N_B}{B}_i}\end{array}\right.$$

(18)

In Formulae (18), $N_A$ and $N_B$ are the number of point clouds in A and B, respectively. The number of points in two point clouds is required to be equal, as $N_A=N_B$.

(2) Calculation the covariance matrix of A and B:

$${\displaystyle R\left(x\right)=\frac{1}{n}\sum _{i=1}^n\left[(A_i-Z_A){(B_i-Z_B)}^T\right]}$$

(19)

The covariance matrix obtained by Formula (19) can form a 4 × 4 symmetric matrix, $Q\left(R\right)$:

$$Q\left(R\right)=\left[\begin{array}{cc}tr\left(R\right(x\left)\right)& {\Delta }^T\\ \Delta & R\left(x\right)+R{\left(x\right)}^T+tr\left(R\left(x\right)\right)I_3\end{array}\right]$$

(20)

where $tr\left(R\right(x\left)\right)$ is the trace of covariance matrix $tr\left(R\right(x\left)\right)$; $\Delta ={\left[A_{23}{A}_{31}{A}_{12}\right]}^T$, and $A_{ij}={(R\left(x\right)-R{\left(x\right)}^T)}_{ij}$; $I_3$ is a unit matrix of 3 × 3.

(3) The eigenvalues and eigenvectors of the symmetric matrix $Q\left(r\right)$ are calculated and the eigenvectors corresponding to the largest eigenvalues are the rotation transformation vector $q_R={\left[q_0{q}_1{q}_2{q}_3\right]}^T$. It should be noted that in the unit quaternion, $q_0>0$, and $q_0^2+q_1^2+q_2^2+q_3^2=0$. The rotation matrix $R\left(q_R\right)$ can be calculated from $q_R$:

$$R\left(q_R\right)=\left[\begin{array}{ccc}{q}_0^2+q_1^2-q_2^2-q_3^2& 2(q_1{q}_2-q_0{q}_3)& 2(q_1{q}_3+q_0{q}_2)\\ 2(q_1{q}_2+q_0{q}_3)& q_0^2-q_1^2+q_2^2-q_3^2& 2(q_2{q}_3-q_0{q}_1)\\ 2(q_1{q}_3-q_0{q}_2)& 2(q_2{q}_3+q_0{q}_1)& q_0^2-q_1^2-q_2^2+q_3^2\end{array}\right]$$

(21)

(4) After obtaining the rotation matrix, the translation matrix $q_T$ can be solved:

$$q_T=Z_A-R\left(q_R\right)Z_B$$

(22)

The rotation matrix $R\left(q_R\right)$ and the translation matrix $q_T$ are the transformation relationship between the source ski-resort datasets and the target ski-resort datasets. The rotation and translation matrix can be used to adjust the three-dimensional coordinates of the source point cloud, which provides a better initial position for the subsequent fine registration.

3.1.5. Overall Process of Coarse Registration Algorithm

Step 1: Identify the common overlaps A and B between the ski-resort point cloud datasets.

Step 2: The common points of A and B are extracted by the SIFT algorithm.

Step 3: Calculate the normal vectors of the feature points and unify their direction.

Step 4: Calculate the Euclidean distance between the normal vector of a feature point in A and the normal vector of all feature points in B, and select the two closest points. If the distance ratio is less than the threshold, the two points with the smallest distance will be matched.

Step 5: According to the spatial topological relationship of the point cloud, the included angles of the normal vectors between the matched feature points are calculated. According to the mathematical relationship, if it is larger than the set threshold, it will be regarded as the mismatching point and be rejected.

Step 6: Use the unit quaternion method to calculate the rotation matrix and translation vector. If the calculated quaternion does not meet $q_0>0$ and $q_0^2+q_1^2+q_2^2+q_3^2=0$, return to step (2) until the requirements are met.

Step 7: Apply the obtained rotation and translation matrix to the source ski-resort point clouds.

3.2. Fine Registration

After coarse registration, the point clouds of the ski-resort datasets have acquired a good initial position. On this basis, ICP iterative algorithm is introduced to improve the effect of fine registration [37,38]. The ICP algorithm calculates the motion parameters (rotation and translation) according to the geometric features between the source ski-resort datasets and the target ski-resort datasets, and then transforms the data with these motion parameters to get the new ski-resort point cloud datasets; in order to continue to determine the registration point cloud between the new correspondence, repeat the above process. The rotation and translation of ski-resort datasets will not stop until the minimum objective function meets the least square theorem. At this time, the registration effect is the best. The objective function is expressed as follows:

$${\displaystyle f(T·R)=\frac{1}{n}\sum _{i=1}^n{\parallel P_i-(Q_i·T·R)\parallel }^2}$$

(23)

where, T and R represent translation and rotation parameters, respectively; $P_i$ and $Q_i$ represent the target ski-resort datasets and source ski-resort datasets, respectively.

4. Results and Discussion

In order to verify the effectiveness of the above algorithm for the point cloud registration of ski resort, this paper takes two groups of laser-point cloud datasets from different locations of the Wanlong ski resort for the registration test. The algorithm in this paper is implemented by C++, and the program running environment is Intel core i5-9400f, DDR4 16 GB, 1080Ti and Windows 10. In order to quantitatively evaluate the registration accuracy of the algorithm, the root mean square error (RMSE) is defined to represent the difference between the features of the same name. The calculation formula is as follows:

$$R_e=\sqrt{\frac{1}{N}\sum {\parallel p_i-q_i\parallel }^2}$$

(24)

In Formula (24), N is the number of overlapping points of two points to be registered after Euclidean transformation; $p_i,q_i$ are the points of the overlapping region after Euclidean transformation.

In the verification process, the two groups of ski resort point clouds datasets are registered, respectively, and the threshold of feature point matching is set to 0.05, and the threshold of matching point pair purification is set to 0.8. The method proposed in this paper and ICP are used to register points, respectively; the registration results are shown in Figure 4 and the rotation matrix and translation vector are shown in

Table 1. The rotation and translation matrix of the first group of ski-resort point clouds calculated based on the method proposed in this paper.

	Rotation Matrix				Translation Vector
R	0.999	−0.010	0.005	T	4.644
	0.010	0.999	0.002		−4.868
	−0.005	−0.002	0.999		3.819
	0	0	0		1.000

.

According to the intersity information in the scalar field of the snow point clouds and the specific color scale, the experimental point cloud datasets are colored into different colors. Figure 4 represents the position distribution among the source point clouds in the first group of Wanlong ski resort, and Figure 5a is the registration results using the ICP, where red marked areas indicate poor registration results with significant obliquity and deviation, Figure 5b is the registration result by the method proposed in the paper. Additionally, the rotation matrix R and translation matrix T are also calculated, the Table 1 is the result of R and T, using the method proposed in this paper.

Figure 6 is the position distribution among the source point clouds in the second group of Wanlong ski resort, Figure 7a is the registration results using the ICP, where blue marked areas indicate poor registration results with significant obliquity and deviation; Figure 7b is the registration result by the proposed algorithm.

In addition, the rotation matrix R and translation matrix T are also calculated; the

Table 2. The rotation and translation matrix of the second group of ski-resort point clouds calculated based on the method proposed in this paper.

	Rotation Matrix				Translation Vector
R	1.010	−0.010	0.007	T	119.057
	0.010	1.010	0.003		78.767
	−0.007	−0.003	1.010		−3.616
	0	0	0		1.000

is the result of R and T using the method proposed in this paper. The actual registration results demonstrate that this group of rotation translation matrix can realize the effective registration between the test ski-resort point cloud datasets.

Then, the registration results of the proposed method and ICP method are compared in terms of the number of point cloud datasets, processing time and RSME. The results are shown in the

Table 3. Comparison of the registration results between the proposed method and ICP method.

Name of Points	Number of Points	Process Time /s		RSME /m
		ICP Method	Proposed Method	ICP Method	Proposed Method
Group one	234,644	61.526	16.951	1.619	0.511
Group two	212,977	41.235	10.962	1.596	0.447

. It can be observed from the Table 3 that the number of point clouds of the first group of ski resort is 234,644 points. Under the same hardware conditions, the processing time of ICP method is 61.526 s, while the method proposed in this paper only takes 16.951 s; similarly, the RMSE of this group of data processed by ICP method is 1.619 m, while the RMSE obtained by our method is only 0.511.

In Figure 6 and Figure 7, when comparing the registration results of the ICP method and proposed method, we can find that there is obvious deviation in the results of ICP registration, and the results of the normal vector registration is better. This proves that the method proposed in this paper is more accurate than the ICP method in registration. In addition, in Table 1, it can be easily found that when processing the same set of data, the method proposed in this paper is more efficient in both the processing time and mean square error than the direct use of the ICP method. In conclusion, the proposed approach on the registration precision and registration efficiency is superior to directly using the ICP algorithm.

In addition, in order to further verify the effectiveness of the algorithm, the comparison is carried out from the dimension of the real error, that is, the real distance between the marked points is compared horizontally with the distance between the marked points in the registered three-dimensional model, and the measurement error is recorded. Taking the second group of ski-resort point cloud datasets as an example, five marker points are selected, namely A, B, C, D and E, and their distribution in the image and the point cloud is shown in Figure 8.

For the five marked points in Figure 8, two methods of field measurement with the meter ruler and three-dimensional measurement with point cloud processing tool are used successively.

In order to reduce the influence of error factors, the interval sampling measurement is carried out between the marked points (A, B, C, D and E) and the target marked points in turn. The measurement results are counted in Table 2 (here, it is considered that the field measurement value with the meter ruler is the true value).

Table 4. Statistical table for error analysis of measuring distance between mark points in different ways.

Marked Points	Target Points	True Distance /m	Measured Distance /m	Error Value /m
A	C	17.224	17.252	0.028
	E	37.805	37.844	0.039
B	D	18.462	18.456	0.006
	E	30.868	30.835	0.033
C	A	17.224	17.270	0.046
	E	20.725	20.698	0.027
D	A	38.225	38.252	0.027
	B	18.462	18.432	0.030
E	A	37.805	37.839	0.034
	C	20.725	20.701	0.024

shows the error analysis and statistics of measuring the real distance, and measured the distance between marked points by tape measurement and point cloud measurement. From the analysis of the results in Table 4, it can be observed that with the increase in measuring distance, the proportion of measuring error will be reduced to a certain extent, and the error range is basically within 5 cm; the good adaptability and anti-noise performance of the proposed method in dealing with the point cloud in ski resorts are further verified.

5. Conclusions and Future Work

According to the complex environment of the Wanlong ski resort, this paper sets the threshold of the Euclidean distance and included angle between the adjacent normal vectors of point clouds, which can effectively accomplish the valid registration of the cloud point datasets of the ski resort. This composite registration method can extract the feature points by SIFT, pair them according to the Euclidean distance between normal vectors of points, set the included angle threshold for further purifying the feature point pairs, calculate the rotation and translation matrix, and finally, complete the fine registration with ICP. The results demonstrate that the proposed composite registration algorithm on paper is more accurate and effective than the direct ICP. The future work will focus on the research of the fine registration method of the point cloud datasets of ski resort, and expect to complete the automation of the whole registration process. Additionally, we will further study the snow characteristics to improve the theoretical foundation for the introduction of more efficient registration methods.