AUTOMATIC GLOBAL REGISTRATION BETWEEN AIRBORNE LIDAR DATA AND REMOTE SENSING IMAGE BASED ON STRAIGHT LINE FEATURES

An automatic global registration approach for point clouds and remote sensing image based on straight line features is proposed which is insensitive to rotational and scale transformation. First, the building ridge lines and contour lines in point clouds are automatically detected as registration primitives by integrating region growth and topology identification. Second, the collinear condition equation is selected as registration transformation function which is based on rotation matrix described by unit quaternion. The similarity measure is established according to the distance between the corresponding straight line features from point clouds and the image in the same reference coordinate system. Finally, an iterative Hough transform is adopted to simultaneously estimate the parameters and obtain correspondence between registration primitives. Experimental results prove the proposed method is valid and the spectral information is useful for the following classification processing. * Corresponding author lpclqq@163.com


INTRODUCTION
Since LiDAR data and remote sensing image are obtained by various sensors under different environmental conditions, these two data have distinct imaging mechanisms, data properties and geographical references.The accuracy of registration between LiDAR and images determines the accuracy of the LiDAR data subsequent processing.In addition, the high-precision registration between LiDAR data and images can simultaneously obtain the spatial and semantic information of the objects, and effectively make up for the shortcomings of a single data source (Dowman, 2004;Baltsavias, 1999).And it is widely used in many applications such as orthophoto generation, disaster assessment, road extraction and building modeling (Smith, et al., 2004).
The registration between LiDAR data and images is part of the registration between continuous image and discrete point could, two-dimensional and three-dimensional data.Traditional registration methods are not suitable for it.In recent years, scholars at home and abroad have done a certain amount of research on them, which can be roughly classified into three types of algorithms:  The corresponding relationship is established betweeen LiDAR data and the three-dimensional point cloud in the object space, which are genenrated by matching multiple aerial images in the same viewpoint via image mathing method.Such methods are essentially a 3D-3D registration transformation.Typically, such methods require better initial values and they can not be used for the registration between LiDAR data and one single image.In addition, the accuracy of point cloud, which is generated via dense matching of images, is generally lower, which reduces the registration accuracy.
 The intensity image or the distance image is generataed by interpolating the LiDAR data, and the registration between LiDAR distance image or intensity image and remote sensing image is realized via traditional image registration methods.The advantages of the above methods are they can make full use of mature image registration methods and have a high degree of automation.The disadvantage is the eror in the process of interpolation and inaccurate point selection, which can reduce the final registration accuracy.
 These registration methods are implemented by directly using LiDAR data and remote sensing image.The corresponding features are extracted from LiDAR data and remote sensing image.The most common features are linear and planar features.The advantages of such methods are they do not necessitate to grid LiDAR data and avoid the dense matching of remote sensing images.The disadvantage is the automatic realization of such methods depends on the automatic extraction and matching of registration primitives, and now the degree of automation needs further research.
In view of the advantages and disadvantages of these abovementioned algorithms, an automatic global registration between airborne LiDAR data and remote sensing image based on straight line features is proposed.

The Automatic Extraction of LiDAR Line Features based on the Fusion of Region Growth and Topology Recognition
The automatic extraction of LiDAR line features is proposed which is based on the fusion of region growth and topology recognition.Firstly, the LiDAR linear features are automatically detected as registration primitives by integrating region growth and topology identification.The point clouds of building area are automatically detected by combining the height threshold and normal line.The point clouds belonging to distinct buildings are separated via the euclidean clustering and the roof point cloud is partitioned by regional growth.The building ridge lines are automatically detected by combining Alpha Shaple and RANSAC algorithm, and refined processed by micro-element feature method.According to the topological relationship between the roof spaces, the public cross-line features of the building are analysised in order to obtain the building contour lines.

The Detection of Point Cloud in Building Roof:
In order to detect the point cloud in the roof, the point cloud of the building area need to be extracted firstly.The DEM is obtained via filtering, and the nDSM is gained by subtracting the DEM from the DSM.The nDSM includes buildings, trees, low shrubs, vehicles and other objects, and the topological relations of plane position between some buildings, buildings and trees may be adjacent, so the extraction of building point cloud may be difficult.Firstly, the height threshold (e.g. 3 meters) is applied to nDSM, and non-ground points such as low vehicles and shrubs are filtered out.The remaining laser footpoints are mainly in the area of building and high vegetation.According to the local spatial distribution of footpoints, most vegetation footpoints can be eliminated effectively and the point cloud of building can be obtained.
After obtaining the point cloud data in the building area, the Euclidean distance clustering method is used to divide the point cloud data of different buildings.The clustering refers to the process of dividing a set of objects into multiple subsets composed of similar objects.The common classification statistics include distance and similarity measure.Through the Euclidean distance clustering segmentation method, the three dimensional point cloud data in different regions can be divided into smaller parts, thus greatly reducing the overall processing time.
The point cloud on the roof of the building is deteted based on regional growth.Firstly, the curvature values of all feet are calculated, and all feet are sorted according to the curvature values.Since the feet with minimum curvature values are located in flat areas, the region growing algorithm selects the foot points with minimum curvature values to grow.

2.1.2
The Extraction of the Straight Lines on the Roof of the Building: Since the shape of the roof is complicated, more flexible and practical method is selected for the plane contour extraction of building roof in the lack of prior knowledge.The Alpha Shape algorithm is adopted because it has the ability to effectively extract the edges of the discrete point set.The principle is shown in figure 1 (Wei, et al., 2004).When the Alpha Shape is used to deal with the plane area, the threedimensional discrete point cloud needs to be projected to its two-dimensional plane to facilitate the subsequent processing.
The extracted roof contour usually consists of several straight lines.For any straight line, other straight lines are believed to be noise.The RANSAC algorithm is used to detect the roof contour of the building.
When the spatial resolution of point cloud is low, the building contour does not fall on the roof edge, so the extracted contour is not strictly in the actual edges of the building.In order to improve the extraction precision of the contour, micro-element feature method are used to refined process the contour (Liang, et al., 2014), and the principle of the method is shown in Figure 2. When the absolute value of the scan angle is larger than the threshold, the laser pulse hit the wall of the building and forms the point cloud.The roof contour of the building and its projection on the ground can be composed of vertical surface.
The element features mean the micro variable of vertical surface at this characteristic.When the number of footpoint within the element features satisfies the threshold, the wall facade can be fitted according to the least square adjustment, and then the outline of the roof contour is refined based on the horizontal position of the wall facade.
When the number of point cloud data on the wall facede is few, it is impossible to use the micro feature to process the contour accurately.In subsequent registration process, the unrefinely processed contours are weighted by smaller weights.

2.1.3
The extractin of the ridge line on the building: The extraction of the ridge line is based on the topological relations between every planes of the roof.Firstly, the adjacency relationships between the roofs should be detertemind, and the theoretical basis is the distance between a plane boundary and the rest of the plane boundaries.The method is as follows: Step 1: Any plane i P should be chosen and the middle point of each contour in this plane can be calculated.
Step 2: Another plane ( ) j P j i  is selected, and the distances are calculated between the center point of each contour and all the footpoint in the plane.
Step 3: If any distance satisfies the threshold (less than the three times of average point cloud spacing), then the plane j P is judged as the adjacency surface of i P and marked.
Step 4: The plane j P is traversed and judegd by step 2 to step 3.
Step 5: The plane i P is traversed and its adjacent surface is judged by step 1 to step 4.
After obtaining the judgment of the topological relations between each plane, the extraction of the ridge line can be realized by calculating the intersection line between the adjacent surfaces.

The Automatic extraction of straight line features of remote sensing image
The method is uesd to extract the straight lines of image (Zhiqing, 2012).Firstly, the edges are detected by Canny operator and the edges are connected based on the edge tracking.
In order to facilitate the follow-up registration, it is necessary to select and optimize the extracted straight line features according to the length, direction and distance of the line segment.
The particular cases are as follows: (1) The straight lines are parallel.If multiple straight line features are approximately parallel, and the distance between these lines is very short, then multiple parallel lines are combined.As shown in Figure 3

The collinear condition equation based on the rotation matrix described by unit quaternion
The geometric relationship between the LiDAR point cloud data and the remote sensing image is very complex, and there is a significant transformation in the scale and rotation.Considering the collinear condition equation is established according to the ideal collinear relationship between the projection center and the image point and the ground point, the LiDAR point cloud data can be regarded as the ground point, and the collinear condition equation is chosen as the registration transformation function to describe the mapping relationship between the ground and the image (Baoming, et al., 2008).
The good initial value is necessary when the collinear condition equation is used to solve the registration transformation parameter, and the iteration may not converge if the initial value is bad.In order to enhance the stability of the method and improve the accuracy of subsequent registration, the unit quaternion is used to describe the rotation matrix instead of the Euler angle.The quaternion has no dependence on the initial value, and has strong reliability and stability (Jun, 2007).A quaternion can be composed of 3 imaginary parts, and its representation is as follows: Where, , ij and k are the imaginary unit.
The operation of quaternion and the rotation matrix described by quaternion are discussed in detail (Jun, 2007).The coordinate rotation transformation matrix M , which is composed of unit quaternion, is shown as follows: q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q The exterior orientation azimuth of the image can be calculated directly by (2), which is shown as follows: = the element vaule of 33  M at column j and row i .
When using the quaternion to pose the image pose, the exterior orientation parameters of the image change to X Y Z q q q q .By using photogrammetry theory, the collinear condition equation based on the quaternion can be obtained.
Where ( , , ) X Y Z = 3D coordinates of LiDAR point cloud data ( , ) xy = coordinates of the corresponding points in the image plane coordinate system ( , , The straight line ab is represented by polar coordinate is shown as follows (Tommaselli, et al., 1998): Where ( , ) = the coordinate any point in the image space coordinate system.
The straight line feature CD of LiDAR data is transformed into straight line feature cd on the image plane through collinear condition equation.If CD and ab are corresponding, then ab and cd should be on the same straight line.So through the transformation of collinear condition equation, the endpoint of straight line feature of LiDAR data should be transformed into the image plane coordinate system and the distance between this endpoint and corresponding line feature in the image is zero, which is defined as the similarity measure.
As shown in Figure 4, the coordinate value of the endpoint of CD is ( , , ) X Y Z , which is transformed into the point c on the plane coordinate system.The coordinate value of c is ( ', ') xy and equation ( 8) is as follows: Similarly, the similar equatin can be obtained through the transformation of the other endpoint D .And so on, 2n equations can be obtained through n pairs of matching groups.
So 3 or more pairs of matching groups can get more than 7 equations, which are more than 6 error equations and one constraint condition equation, and the indirect adjustment method with restricted conditions are used to calculate (Lifen, et al., 2010).

THE OPTIMIZATION OF REGISTRATION TRANSFORMATION FUNCTION PARAMETERS
There is large rotation and scale transformation between LiDAR point cloud and the image.It is necessary to find suitable methods for parameter optimization in order to reduce search space and computation complexity.
There are some problems to affect the accuracy and efficiency of parameter calculation when using Hough transform.The main problems are as follows:  The problem of efficiency and memory.Meantime in order to calculate seven exterior orientation parameters, four pairs of matching groups are needed to get eight equations.However, this method will cause the expansion of the number of groups and sharp increase in the amount of calculation; in addition, the 7 dimensional matrix is needed to store data, which also has a great influence on memory.
 The problem of peak search.When searching for peak, the choice of parameter step is very important.If the step size is chosen smaller, it will affect the memory; if the step size is too large, it will reduce the accuracy of peak selection.
 The elimination of false matching groups.When the error equation is established in order to calculate the exterior orientation parameters by using the matching groups contributing to the peak, the false matching groups will have a great impact on the accuracy of the parameters and registration.
Aiming at the above-mentioned problems, some solutions are proposed including the reduction the matrix dimension, pyramid search, the iteration method with variable weights and the judgment of distance.
 The reduction of the matrix dimension.A cumulative matrix is established for 7 parameters respectively.The parameters are calculated by making use of assumed matching groups, and the peak value is obtained by accumulating the corresponding accumulation matrix.Finally, the accumulative matrix also becomes a one-dimensional array, and the memory problem is solved.
 The pyramid search.The pyramid principle is used to solve the problem of setting the peak search step.For the first time, the step size should be larger in order to avoid the elimination of adjacent peak point.Then, the calculated parameters are stored in the accumulator array in order for frequency statistics; the interval with maximum frequency is expanded 1/2 step on the left and right sides, which is set to the new range of accumulator array.If the size of new range meets the threshold, then the iteration stoped; or the iteration contimues in order to obtain the new range until the new range meets the threshold.
 The iteration method with variable weights and the judgement of distance.Firstly, the Danish method is introuduced to the calculation of exterior orientation parameters to eliminate the errors with the iteration method with variable weights, and so as to improve the accuracy and reliability of parameter calculation.The judgment of distance is used to further eliminate the false matching groups.The matching groups contributing to the peak are used to calculate the parameters, and the straight line features of LiDAR from matching groups are projected to the image via the calculated parameters.Through the factors such as the direction, distance, and whether overlapping, the false matching groups are eliminated.A specfic example is shown in Figure 5.The straight line ab is on the image, and respectively with straight lines 12, 34, and 56 of LiDAR data ab are matching groups contributing to the peak.These straight lines 12, 34, and 56 are projected on the image and firstly the straight line 56 is judged since the angle is larger between ab and 56.The straight lines ab is approximately parallel with 12 and 34 respectively, and the distance between 34 and ab is shorter than that of 12 and ab .There is an overlap between 12 and ab , but no overlap between 34 and ab , and whether whether there is some overlap is the most important criteria.So the straight line 34 is eliminated and 12 is reserved.
The flow of the optimal method is as follows: Step 1: The pair of matching group is randomly selected, and the angle between the lines of group is calculated.If the angle is less than the threshold, the next step contimues, otherwise new groups are chosen.
Step 2: The error equation is established to calculate the parameter s X , and according to the calculated value the corresponding interval is to vote for.
Step 3: The step 1 and 2 are repeated until all hypothesis matching combinations are traversed.
Step 4: The peak value is selected for parameter s X , and the mean of peak value is as the initial value of parameter for the next iteration.The matching groups contributing to the peak value are recorded.
Step 5: The other six position and attitude parameters are calculated according to the step 1 to 4 in the same way, and when seven parameters are all calculated we can start step 6.
Step 6: The new initial values are for seven parameters and the contributing matching groups are computed.The error equation is established by using these matching groups and the calculation accuracy is improved through the iteration with variable weights.If the residual is less than the threshold, the next processing starts; or the initial value is set as the parameter value for last iteration, and the recalculation starts from step 1.
Step 7: The matching groups processed by step 6 are filtered again according to the factors such as direction and distance, and the matching groups after screening are used to establish the error equations so as to calcuate the parameters.

EXPERIMENT AND ANALYSIS
In order to prove the performance of the automatic registration method between LiDAR data and remote sensing image based on straight line features, the data set from Vaihingen and African area were selected for experiments and the experimental results are illustrated and analyzed finally.From the point cloud rendering results, some conclusions can be found: (1) the proposed registration method is good; (2) combined with spectral information, the characteristics of discrete point cloud objects are more obvious and easy to identify, indicating that spectral information will provide important help for identifying objects.

The experiment and analysis of data in Vaihingen
Besides the qualitative experiments, the quantitative experiments are also needed to verify the registration accuracy, which includes checking residual error statistics.The seven straight line features are extracted from LiDAR data as the checking lines, which are projected on the image.The measurement standards are the mean square error of the maximum distance, the minimum distance, the average distance and the angle error respectively between the projected LiDAR line features and the corresponding line features of the image.The specific quantitative experimental results are shown in table 1.The mean square error of the maximum distance, the minimum distance, the average distance and the angle error are 1.84 pixels, 1.31 pixels, 1.57 pixels and 1.35 degrees respectively between the projected LiDAR line features and the corresponding line features of the image.Since the resolution of the image is 0.08 meters, the mean square error of the maximum distance, the minimum distance and the average distance can be converted into 0.15 meters, 0.11 meters and 0.13 meters respectively.The point spacing of LiDAR data is about 0.4 meters, and the elevation accuracy is about 0.15 meters.The specific quantitative experimental results are shown in table 2. The mean square error of the maximum distance, the minimum distance, the average distance and the angle error are 1.66 pixels, 0.71 pixels, 1.15 pixels and 1.14 degrees respectively.Since the resolution of the image is 0.1 meters, the mean square error of the maximum distance, the minimum distance and the average distance can be converted into 0.17 meters, 0.07 meters and 0.12 meters respectively.The point spacing of LiDAR data is about 0.55 meters, and the elevation accuracy is about 0.15 meters.The qualitative and quantitative experimental results meet the accuracy requirements. No

CONCLUSION
In this paper, an automatic global registration between airborne LiDAR point cloud data and remote sensing images based on straight line features is researched.This method is insensitive to rotation and scale changes between point cloud data and remote sensing images.The experimental results show that the proposed method is valid, and the spectral information provides an important reference for the subsequent classification.

Figure 1 .
Figure 1.The priciple of Alpha Shape algorithm (a), AB and CD are combined into A'B'.(2) The straight lines are overlapped.If multiple straight linear features are overlapped, only the longest straight line is reserved.As shown in Figure 3 (b), among the overlapped AB, CD and EF, AB is only reserved.(3) The straight lines are fracture.If one straight line is divided into several small segments, the longest straight line segment is reserved.As shown in Figure 3 (c), among the broken AB, CD and EF, AB is only reserved.The selection and optimization of line features

Figure 4 .
Figure 4.The similarity measure According to the above-mentioned figure, through the conversion based on the collinear condition equation, the straight line feature CD of LiDAR data becomes the straight line feature cd in the image plane coordinate system, and cd

Figure 5 .
Figure 5.The principle of the judgment of distance

Firstly
, the straight lines of LiDAR are extracted automatically.The specific case is as follows: (a) is the automatic detection of the building; (b) is the point cloud clusters of different buildings; (c) is the roof detection of the building; and (d) is the ridge lines and contour of the building.(a) the point cloud detection of building (b) the point cloud clustering of building (c) the point cloud detection on the roof (d) the extraction of the ridge lines and contour lines Figure 6.The straight line features extraction of LiDAR point cloud in VaihingenThe optimized extraction of line features on the remote sensing image is shown in figure7.The point cloud of single building and tree are selected and projected onto the image through calculated registration transformation parameters.The specfic case is shown in figure8.It can be seen that the registration results have no obvious dislocation and meet the requirements of subsequent classification processing.

Figure 7 .
Figure 7.The optimized extraction of line features on the remote sensing image

Figure 9 .
Figure 9.The point cloud rendering results conbining with spectral information in Vaihingen

5. 2
The experiment and analysis of data in African area The straight line features extraction of LiDAR data in African area is shown in figure 10.(a) the point cloud detection of building (b) the point cloud clustering of building (c) the point cloud detection on the roof (d) the extraction of the ridge lines and contour lines Figure 10.The straight line features extraction of LiDAR point cloud in African area The optimized extraction of line features on the remote sensing image is shown in figure 11.The point cloud of single building and tree are selected and projected onto the image, which are shown in figure12.It can be seen that the registration results have no obvious dislocation and meet the requirements of subsequent classification processing.

Figure 11 .
Figure 11.The optimized extraction of line features on the remote sensing image

Figure 13 .
Figure 13.The point cloud rendering results conbining with spectral information in African area

Table 1 .
The checking line residual in Vaihingen The qualitative and quantitative experimental results show the proposed method is valid.

Table 2 .
The checking line residual in African area