Registration of Urban Aerial Image and LiDAR Based on Line Vectors

In a traditional registration of a single aerial image with airborne light detection and ranging (LiDAR) data using linear features that regard line direction as a control or linear features as constraints in the solution, lacking the constraint of linear position leads to the error propagation of the adjustment model. To solve this problem, this paper presents a line vector-based registration mode (LVR) in which image rays and LiDAR lines are expressed by a line vector that integrates the line direction and the line position. A registration equation of line vector is set up by coplanar imaging rays and corresponding control lines. Three types of datasets consisting of synthetic, theInternational Society for Photogrammetry and Remote Sensing (ISPRS) test project, and real aerial data are used. A group of progressive experiments is undertaken to evaluate the robustness of the LVR. Experimental results demonstrate that the integrated line direction and the line position contributes a great deal to the theoretical and real accuracies of the unknowns, as well as the stability of the adjustment model. This paper provides a new suggestion that, for a single image and LiDAR data, registration in urban areas can be accomplished by accommodating rich line features.


Introduction
LiDAR is characterized by high elevation accuracy, strong autonomy, and well-operated automatization [1], which could directly obtain the three-dimensional information of ground features.The discreteness in LiDAR data makes it difficult to directly obtain semantic information (e.g., texture and structure) of ground features [2].Aerial images are abundant in spectral information and texture information of ground features, therefore, the surface accuracy is much higher.However, there are potential risks in three-dimensional points with sufficient density acquired by auto-matching superimposed images.LiDAR data and aerial images are highly complementary [2], so the combination of these two sources of data could simultaneously obtain three-dimensional, semantic, and texture information from spatial targets.This has important significance for feature extraction [3], three-dimensional reconstruction of buildings [4][5][6], and manufacturing of true orthophotos [7].Various reference frames in airborne LiDAR data and imagery, along with systematic errors inside of the LiDAR system and aerial photography system, cause spatial position deviations between the LiDAR data and the image, which may greatly affect the accuracy of the data results.Both of them, before being integrated for application, should be included in a unified coordinate system to achieve their spatially-accurate registration, namely the registration of LiDAR data and aerial images [8].
The commonly-used features in registration mainly include point features and the line features.Compared with point features, linear features have the potential to be reliably, accurately, and automatically extracted from images [9].In addition, linear features are abundant and can be easily obtained in urban areas without the requirement for the line segment endpoint corresponding to the same point, allowing the linear features to be partially occluded [10][11][12][13][14][15].
In registration based on linear features, Habib et al. proposed to directly incorporate LiDAR lines representing a sequence of intermediate points along the line features as controls in bundle adjustment [16,17].Liu et al. proposed to establish coplanar equations of bundle adjustments by the coplanar constraints of linear features and collinear constraints of tie points [18].However, this requires appropriate initial exterior orientation parameters (EOPs) to obtain approximate locations when projecting three-dimensional (3D) line segments in the LiDAR data coordinate framework to the image coordinate framework.These methods solve the position parameters and direction parameters separately when solving EOPs due to which coupling error exists.To overcome the probable error accumulation in the process of asynchronous calculation of transformation parameters, dual quaternions are used to describe the spatial rotation and to ensure the simultaneous calculation of registration parameters.However, the registration model of this method is a spatial similarity transformation, whose registration primitive is a feature point [19,20].
Almost all of the aforementioned methods use the collinearity equation or coplanarity condition as a mathematical foundation for registration and regard line direction as a control or linear features as known variables in the solution, whose adjustment form is identical to the point registration method.In fact, any line in the space is solely determined by its own position and direction [21,22].If merely using the direction information for registration, the geometry constraint of registration model is not strong enough.
On the basis of the existing methods with linear features, registration of an urban aerial image and LiDAR based on line vectors that aggregate line direction and line position is proposed in this paper.The remainder of this paper is organized as follows: The principle of line vectors is first introduced.Then, the mathematical model of the registration method based on line vectors is discussed in detail.Finally, we present the experimental results and our conclusion.

Methods
The registration based on line vectors involves two steps.The first step is, respectively performing line description and other operations.The second step is to apply the two sets of data obtained in the first step to the coplanarity equation and find the EOPs satisfying the iteration condition.The overall flow chart in this paper of registration based on line vectors is shown in Figure 1.

Line Vector
As shown in Figure 2, the unit direction vector of lines can be expressed by three direction cosines, namely = ( , , ) , which satisfies + + = 1 , i.e., = 1 , but these three parameters cannot determine the position of the line.The position vector of H, an arbitrary point on the line directed by the origin is r, according to which the position of the line can be determined and line vector L can be obtained.Line vector L is a six-vector, including the primary part and secondary part.The directional vector is the primary part of the line vector, while the moment vector is the secondary part of the line vector.The directional vector does not change as the origin changes, but the moment vector does.For this reason, we know that the line vector is a vector attached to a line in the space and six-dimensional linear space in the real number field: The directional vector and moment vector satisfy: According to the two constraints that = 1 and • = 0, which derive from the unit line vector itself, there are four degrees of freedom in the line vector coordinate.

Line Vector
As shown in Figure 2, the unit direction vector of lines can be expressed by three direction cosines, namely o = (l, m, n) T , which satisfies l 2 + m 2 + n 2 = 1, i.e., o T o = 1, but these three parameters cannot determine the position of the line.The position vector of H, an arbitrary point on the line directed by the origin is r, according to which the position of the line can be determined and line vector L can be obtained.Line vector L is a six-vector, including the primary part and secondary part.The directional vector is the primary part of the line vector, while the moment vector is the secondary part of the line vector.The directional vector does not change as the origin changes, but the moment vector does.For this reason, we know that the line vector is a vector attached to a line in the space and six-dimensional linear space in the real number field: The directional vector and moment vector satisfy: According to the two constraints that o T o = 1 and p • o = 0, which derive from the unit line vector itself, there are four degrees of freedom in the line vector coordinate.) .Operations on line vectors are performed similar to those on Cartesian vector, also referred to as line dot-and cross-products [17], and are defined by: The advantage of the linear expression method based on line vector coordinate is in determining a line not relying on any point or surface, as well as directly representing a line direction by using a directional vector and line position by a moment vector.
Assuming any two lines with their directional vectors set as and , respectively, and their moment vectors as and , respectively, if these two lines are coplanar, then they will meet the requirement of: Equation ( 4) is the coplanar formula of the two lines in a simple way.Compared with the two traditional methods of using point coordinates and the orientation vector to express the coplanar formula, the line vector may directly make use of the information of direction and position.

The Coordinate Transformation of Line Vectors
As shown in Figure 3, hg is an image feature line passing through point h and point g.Imaging ray Sh from the photographing center S runs through point h and Sg runs through point g.As a corresponding control line in LiDAR data, GH runs through point G and H.In regard to lines Sh, Sg, and GH, assume their directional vectors of are , , and O, moment vectors are , , and P, and all three lines are coplanar at photographing, then, based on Equation ( 4), the LiDAR data coordinate system and the image coordinate system have the following relation: In the aforesaid equation, the coordinate transformation relationship between the LiDAR data coordinate system and the image coordinate system is characterized by representing coordinates of g and h in the image coordinate system, as well as coordinates of G and H in the LiDAR coordinate system.Any two LiDAR points can determine control lines and two image points can determine corresponding image lines, so g and h, as well as G and H, are not required to be corresponding points, which is the advantage of registration using line characteristics.Assume there are two line vectors L 1 and L 2 , which are Operations on line vectors are performed similar to those on Cartesian vector, also referred to as line dot-and cross-products [17], and are defined by: The advantage of the linear expression method based on line vector coordinate is in determining a line not relying on any point or surface, as well as directly representing a line direction by using a directional vector and line position by a moment vector.
Assuming any two lines with their directional vectors set as o 1 and o 2 , respectively, and their moment vectors as p 1 and p 2 , respectively, if these two lines are coplanar, then they will meet the requirement of: Equation ( 4) is the coplanar formula of the two lines in a simple way.Compared with the two traditional methods of using point coordinates and the orientation vector to express the coplanar formula, the line vector may directly make use of the information of direction and position.

The Coordinate Transformation of Line Vectors
As shown in Figure 3, hg is an image feature line passing through point h and point g.Imaging ray Sh from the photographing center S runs through point h and Sg runs through point g.As a corresponding control line in LiDAR data, GH runs through point G and H.In regard to lines Sh, Sg, and GH, assume their directional vectors of are o h , o g , and O, moment vectors are p h , p g , and P, and all three lines are coplanar at photographing, then, based on Equation ( 4), the LiDAR data coordinate system and the image coordinate system have the following relation: In the aforesaid equation, the coordinate transformation relationship between the LiDAR data coordinate system and the image coordinate system is characterized by representing coordinates of g and h in the image coordinate system, as well as coordinates of G and H in the LiDAR coordinate system.Any two LiDAR points can determine control lines and two image points can determine corresponding image lines, so g and h, as well as G and H, are not required to be corresponding points, which is the advantage of registration using line characteristics.Directional vectors and moment vectors of imaging rays cannot make Equation ( 5) workable because there are certain deviations of the original position and direction of the photographing center.Such a photographing center should be moved to restore the correct position and direction of Sh and Sg.That is to say, Sh and Sg should be carried out of the coordinate transformation to meet Equation ( 5).The coordinate transformation of the line vector is generally achieved by a dual quaternion.Line can be obtained from L by the dual quaternion, q: = where and is the conjugation of q.

Coplanarity Equation Based on Line Vectors
Assuming the original line vector coordinate of the imaging ray Sh is (l 1 , m 1 , n 1 , r 1 , s 1 , t 1 ).Then the line vector coordinates of Sh through coordinate transformation can be obtained by Equation ( 6): Similarly, line vector coordinates of Sg through coordinate transformation can be obtained by its original line vector coordinate (l2, m2, n2, r2, s2, t2).Assuming the line vector coordinate of GH is (L, M, N, R, S, T), wherein O = (L, M, N) and P = (R, S, T).The coplanarity equation based on line vectors can be expanded by substituting the line vector coordinates of Sh, Sg, and GH into Equation (5).Sh, Sg, and GH are coplanar when Sh and Sg obtain the correct position and direction by correcting Directional vectors and moment vectors of imaging rays cannot make Equation ( 5) workable because there are certain deviations of the original position and direction of the photographing center.Such a photographing center should be moved to restore the correct position and direction of Sh and Sg.That is to say, Sh and Sg should be carried out of the coordinate transformation to meet Equation ( 5).The coordinate transformation of the line vector is generally achieved by a dual quaternion.Line L can be obtained from L by the dual quaternion, q: where q = q 1 + q 2 i + q 3 j + q 4 k + ε(q 01 + q 02 i + q 03 j + q 04 k), i 2 = j 2 = k 2 = ijk = −1, ε 2 = 0, and q −1 is the conjugation of q.

Practical Procedure of Solving the Problem
The proposed adjustment mode includes the following five parts, and all the data are associated automatically: (1) Representation of a Straight Line The unit line vector coordinate constituted by any two points G (X G , Y G , Z G ) and H (X H , Y H , Z H ) of a control line in the LiDAR data coordinate system is [19]: Similarly, line vectors can be used to express two imaging rays Sh and Sg. (

2) Determination of Initial Values
We can determine the original values of the unknown values q by inputting rotation parameters ω, φ, κ, and translation parameters X L , X L , X L into Equation ( 9): (3) Establishment of the Error Equation For the least-squares adjustment, the coplanarity equations (Equation ( 5)) are linearized in the form: where the observations are the image coordinate measurements, and the unknowns are the registration parameters described by X. V are the correction values of the observed values.F are the residual terms of the observed values, which are obtained by putting the approximate value of the unknown value on the left-hand side of Equation (5).A are the coefficient terms of the error equation.

Testing Methods
Three datasets of corresponding lines from the synthetic, ISPRS test project, and real image are adopted for experiments.The synthetic image is adopted to verify the correctness of the proposed algorithm, and the ISPRS test project image and real image are adopted to verify the practicability of the proposed algorithm.Conventional registration and the proposed registration method based on line vectors are performed in each type of image.The initial set of unknowns in conventional registration is X L = Y L = Z L = ω = φ = κ = 0, while the LVR gives data of q 1 = 1, q 2 = q 3 = q 4 = q 01 = q 02 = q 03 = q 04 = 0, equivalent to the initial set in conventional registration.A description of each dataset is described below.
(1) Synthetic images: The generation method of synthetic data is as follows: firstly a parameter distortion model (PDM) is used in the camera calibration with radial distortion and decentering distortion included.Then the accurate camera inner orientation parameters (IOPs) and the EOPs of each image are obtained.Substituting the ground control points, camera parameters, and the EOPs of the image into the collinear equation, the coordinates of corresponding points can then be calculated.Finally, the IOPs and the corresponding points in the image and LiDAR data coordinate framework are used as the given data to solve the EOPs, since there may be some random errors in the measurement of point coordinates in practical applications.To simulate the actual situation, we using the Gaussian random function to produce Gaussian errors with zero-mean and 1.0 pixel standard deviation based on simulated data, assign these data in the x-axis and y-axis of 18 endpoints of image lines, wherein random errors on the x-axis and y-axis are required to be, respectively, generated to ensure their independence.The IOPs of the synthetic image have a principal point of x 0 = 0, y 0 = 0 and a principal distance of f = 28 mm.The image format was 1024 × 1280 with a pixel size of 0.008 mm.
(2) ISPRS test project data: A group of data is a subset of the ISPRS test project on Urban Classification and 3D Building Reconstruction.Everyone can download it from the ISPRS website.The study area is located in Vaihingen City.Aerial images are shot by a digital mapping camera (DMC), wherein the focal length of DMC is 120 mm, the charge-coupled device (CCD) pixel size is 12 µm ground resolution of images and the frame of image is 7680 × 13,824.LiDAR point clouds are collected by the Leica ALS50 subminiature LiDAR with a density of the point clouds of 6.7 point/m 2 .
(3) Real aerial imagery and airborne LiDAR data: The study area is located at the University of Calgary.Aerial imagery is shot by Phase One, wherein the focal length of DMC is 60.6785 mm, the CCD pixel is 6 µm ground resolution of images and the frame of image is 6732 × 8984.The density of the point clouds is 3.5 point/m 2 .
For each dataset, a single image is selected for registration with the LiDAR data.The image lines and LiDAR lines were manually measured.A set of nine evenly-distributed control lines were manually selected for each dataset, while others were used as check lines.Figure 4 shows the distributed control lines for the ISPRS test project, and the real image, respectively.Yellow lines are used to represent control lines and green lines are used to represent check lines.
Appl.Sci.2017, 7, 965 7 of 11 (1) Synthetic images: The generation method of synthetic data is as follows: firstly a parameter distortion model (PDM) is used in the camera calibration with radial distortion and decentering distortion included.Then the accurate camera inner orientation parameters (IOPs) and the EOPs of each image are obtained.Substituting the ground control points, camera parameters, and the EOPs of the image into the collinear equation, the coordinates of corresponding points can then be calculated.Finally, the IOPs and the corresponding points in the image and LiDAR data coordinate framework are used as the given data to solve the EOPs, since there may be some random errors in the measurement of point coordinates in practical applications.To simulate the actual situation, we using the Gaussian random function to produce Gaussian errors with zero-mean and 1.0 pixel standard deviation based on simulated data, assign these data in the x-axis and y-axis of 18 endpoints of image lines, wherein random errors on the x-axis and y-axis are required to be, respectively, generated to ensure their independence.The IOPs of the synthetic image have a principal point of x0 = 0, y0 = 0 and a principal distance of f = 28 mm.The image format was 1024 × 1280 with a pixel size of 0.008 mm.
(2) ISPRS test project data: A group of data is a subset of the ISPRS test project on Urban Classification and 3D Building Reconstruction.Everyone can download it from the ISPRS website.The study area is located in Vaihingen City.Aerial images are shot by a digital mapping camera (DMC), wherein the focal length of DMC is 120 mm, the charge-coupled device (CCD) pixel size is 12μm ground resolution of images and the frame of image is 7680 × 13,824.LiDAR point clouds are collected by the Leica ALS50 subminiature LiDAR with a density of the point clouds of 6.7 point/m 2 .
(3) Real aerial imagery and airborne LiDAR data: The study area is located at the University of Calgary.Aerial imagery is shot by Phase One, wherein the focal length of DMC is 60.6785 mm, the CCD pixel is 6 μm ground resolution of images and the frame of image is 6732 × 8984.The density of the point clouds is 3.5 point/m 2 .
For each dataset, a single image is selected for registration with the LiDAR data.The image lines and LiDAR lines were manually measured.A set of nine evenly-distributed control lines were manually selected for each dataset, while others were used as check lines.Figure 4 shows the distributed control lines for the ISPRS test project, and the real image, respectively.Yellow lines are used to represent control lines and green lines are used to represent check lines.

Experimental Results
The accuracy of registration consists of theoretical accuracy and real accuracy [23].The standard error of unit weight (σ 0 ) of least squares adjustment is used to estimate the theoretical accuracy of the line vector registration.Theoretical accuracy of unknowns is shown in Table 1.To quantitatively verify the real accuracy of line vector registration, we used the root mean square error (RMSE) of the check line, which is calculated by projecting two arbitrary points of each LiDAR check line onto the image using solved EOPs and then calculating the mean distances from projected points to their corresponding image lines.The errors for check lines of conventional registration and LVR methods are listed in Table 2.
Figure 5 represents the effect on the image before and after registration.We extract the point cloud of building roofs in LiDAR data, and project the roof point cloud onto the image by the solved EOPs, where projected roof points are shown in blue.Figure 6 represents the whole and the local 3D effects of the images after registration.Since the LiDAR was pointing vertically down, the points only distribute on the roof of the building, without any distribution on the facade, and the color of the building facade is black after registration.

Discussion
The registration of aerial images and airborne LiDAR data based on linear features are merely

Discussion
The registration of aerial images and airborne LiDAR data based on linear features are merely methods based on points or establishing a registration model using line direction.The coplanarity equation of the line vector is established by coplaning two imaging rays and the corresponding control lines.A registration method using line vectors is raised to solve the traditional weak geometric constraint.
Experimental results show that both conventional registration and LVR can achieve the correct solution.When handing synthetic data, position error of EOPs from both conventional registration and LVR method is within 0.84 m and an attitude error of EOPs is within 0.00466 • .Compared with conventional registration, LVR has 18% higher accuracy in the position parameter and 14% higher accuracy in the attitude parameter.Additionally, the averaged error of check lines by conventional registration is 1.07 pixels, while with the LVR method it is 0.88 pixels.Compared with conventional registration, the LVR method is 18% better in real accuracy.When handling ISPRS test project data, the EOPs are estimated with a σ 0 value of 2.00 pixels via conventional registration, and the RMSE of the check lines at the image scale is 0.66 pixels.In contrast, a σ 0 value of 1.75 pixels is achieved for the LVR, with the RMSE of the check line being 0.54 pixels.When handling real aerial images in Calgary, the LVR is about 20% higher and 24% higher than conventional registration in theoretical accuracy and real accuracy, respectively.
The proposed adjustment is evaluated by a synthetic, ISPRS test project and real data.Compared with the traditional method, the LVR is theoretically and actually more accurate.This is because conventional registration only uses the direction of the line to set up the registration model, whose essence is still the coplanarity of points.The LVR uses directional and moment vectors to describe feature lines, covering both the direction and position information of the lines.
The inaccurate original EOPs before registration has resulted in large deviation in real roof positions and projecting points on the image.The registration of real roof positions and projecting points on an image has been accurately achieved after the registration of aerial images and LiDAR data.It follows that the LVR can achieve an ideal registration accuracy in urban areas with rich linear features.
There are some limitations in the proposed algorithm, though the imaging rays in the image coordinate framework and their corresponding linear control features in LiDAR data coordinate framework can be effectively achieved by line vectors.The raised method achieves accurate extraction of line features by selecting points on aerial images and LiDAR for plan fitting and then plane intersection.Point features selected artificially are assumed on the same plane, so it is actually difficult to ensure that all these points are on the building surface.The contingency may occur in coplanar features as selected points cannot represent the coplanar constraint conditions.

Conclusions
In this paper, we have proposed a new algorithm that uses line vectors to represent imaging rays and control lines.Moreover, it breaks the limitation of the traditional "four-point" coplanar condition by directly describing the "three-line" coplanar condition.The geometrical meaning of the proposed model for the registration becomes more explicit.
According to this study, straight lines are proven to be reliable and effective in registration, and points in the image coordinate framework and LiDAR data coordinate framework on the lines are not necessarily correspondent.To reduce the contingency in extracting line features, further research will be focused on directly using plane features for registration of aerial image and LiDAR data.

Figure 1 .
Figure 1.Steps of the aerial images and LiDAR data registration based on line vectors.

Figure 1 .
Figure 1.Steps of the aerial images and LiDAR data registration based on line vectors.

Figure 2 .
Figure 2. Line vector coordinate of lines.

Figure 3 .
Figure 3. Coplanar model based on line vectors.

Figure 4 .
Figure 4. Distribution of (a) simulated and (b) ISPRS test project.

Figure 4 .
Figure 4. Distribution of (a) simulated and (b) ISPRS test project.

Figure 5 .
Figure 5. Registration effects of the image.

Figure 6 .
Figure 6.3D display of the registration effects.

Figure 6 .
Figure 6.3D display of the registration effects.

Figure 6 .
Figure 6.3D display of the registration effects.

Table 1 .
Registration results of simulated, ISPRS test project, and real images.

Table 2 .
RMSE of check lines of conventional registration and LVR (pixel).