A dual quaternion-based, closed-form pairwise registration algorithm for point clouds

ISPRS Journal of Photogrammetry and Remote Sensing 94 (2014) 63–69

Contents lists available at ScienceDirect

ISPRS Journal of Photogrammetry and Remote Sensing journal homepage: www.elsevier.com/locate/isprsjprs

A dual quaternion-based, closed-form pairwise registration algorithm for point clouds Yongbo Wang ⇑, Yunjia Wang, Kan Wu, Huachao Yang, Hua Zhang Key Laboratory for Land Environment and Disaster Monitoring of SBSM, Jiangsu Key Laboratory of Resources and Environmental Information Engineering, China University of Mining and Technology, Xuzhou, Jiangsu 221008, China

a r t i c l e

i n f o

Article history: Received 30 August 2013 Received in revised form 31 March 2014 Accepted 15 April 2014

Keywords: LiDAR Registration Similarity transformation Rigid transformation Dual quaternions Geodetic datum

a b s t r a c t The representation of similarity transformation in three-dimensional (3D) space, especially of orientation, is a crucial issue in navigation, geodesy, photogrammetry, robot arm manipulation, etc. Considering the large amount of computer resources required by iterative algorithms designed for spatial similarity transformation, the high dependence on initial values of unknown parameters, and the instability of solving transformation parameters for large-angle registration, a closed-form solution for pairwise light detection and ranging (LiDAR) point cloud registration is proposed. In this solution, dual-number quaternions are used to represent the 3D rotation. The relationship between the rotation matrix-based representation of similarity transformation and the dual quaternion-based representation is described first. Considering that the same features from two neighboring stations coincide after pairwise registration, a dual quaternion-based error norm, which is associated with the sum of the position errors, is constructed. Based on theory of least squares and by extreme value analysis of the error norm, detailed derivations of the model and the main formulas are obtained. Once the similarities between the same features from the two neighboring LiDAR stations are constructed, the rotation matrix, the scale parameter, and the translation vector are simultaneously derived. Two experiments are conducted to verify the feasibility and effectiveness of the proposed algorithm. The proposed algorithm has the advantages of simplicity and ease of implementation, making it better than the traditional methods that use matrices to describe spatial rotation. Moreover, it solves the transformation parameters without the initial estimates of unknown parameters, making it better than iterative algorithms. Most importantly, in contrast to unit quaternion-based algorithms, the proposed algorithm solves seven unknown parameters simultaneously. Therefore, it effectively avoids the accumulation of introduced error in calculation and the negative impact from the inappropriate choice of initial values. Ó 2014 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS). Published by Elsevier B.V. All rights reserved.

1. Introduction Reconstructing geographical entities and their neighborhoods efficiently and accurately is a key problem in digital city and is also a core research topic in three-dimensional (3D) geographic information system technology. Since the former US vice president Al Core introduced the concept of digital earth in 1998, the seamless 3D representation of large amounts of available geographical entities in a computer has been a popular topic of research. Among the newly developed instruments designed for the acquisition of location-based data, terrestrial light detection and ranging (LiDAR) ⇑ Corresponding author. Tel.: +86 136 55208823. E-mail addresses: [email protected] (Y. Wang), [email protected] (Y. Wang), [email protected] (K. Wu), [email protected] (H. Yang), [email protected] (H. Zhang).

has been the focus of considerable attention because of its instantaneous and direct acquisition of point cloud. However, given the high spatial complexity of geographical entities and the limited field of LiDAR’s vision, several stations should be established to fully acquire the detailed point cloud of some entities of interest. Therefore, similarity transformation is used to register point clouds from different stations and integrate point clouds. Given the matched points, i.e. pairs of conjugate points, from the two neighbor point clouds as input, many registration approaches have been proposed to calculate the transformation parameters from the nonlinear overdetermined equations in the least-squares sense. Presently, the most frequently used model is the similarity transformation model with seven parameters (i.e., one scale factor, three translation parameters, and three rotation angles). It is also known as the Helmert transformation or the Bursa model, which is used in this study. By using different

http://dx.doi.org/10.1016/j.isprsjprs.2014.04.013 0924-2716/Ó 2014 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS). Published by Elsevier B.V. All rights reserved.

64

Y. Wang et al. / ISPRS Journal of Photogrammetry and Remote Sensing 94 (2014) 63–69

solutions to calculate transformation parameters, available registration approaches can be divided into two categories, namely, analytical methods (Horn, 1987, 1988; Walker et al., 1991; Daniilidis, 1999; Prošková, 2011, 2012) and iterative methods (Besl and Mckay, 1992; Zheng et al., 2008). The main difference between these two categories lies in the requirement of initial estimates of the unknown transformation parameters before registration. Considering that the Bursa model is a nonlinear model, the linearization of multivariate function is necessary for iterative methods. However, it may lead to the loss of accuracy to some extent, and the convergence of the linearized model may be affected by the closeness of the initial estimates and their corresponding true values. Furthermore, the Bursa model is difficult to implement in cases of large rotation angles because the initial estimates are difficult to obtain in some cases, such as absolute orientation in photogrammetry. By contrast, analytical methods may effectively avoid such a case. By using an optimization process, the analytical method directly calculates the similarity transformation parameters. The difference between analytical methods and iterative methods mainly lies in the representation of the 3D rotation. Currently, two operators used in analytical registration methods, namely, orthonormal matrices (Horn, 1988) and quaternions (Horn, 1987), are the two prominent representatives of vector algebra and Clifford algebra. Orthonormal matrices provide a closed-form solution by singular value decomposition (SVD), and quaternions provide a closed-form solution by minimizing a cost function. Compared with matrix-based methods, unit quaternionbased methods use only four parameters, which is less than the six parameters used in matrix-based methods, to represent spatial rotation. Therefore, the demand for computer resources will significantly decrease. This decreased demand has been the focus of considerable attention worldwide. The disadvantages of unit quaternion-based methods lie in the fact that the rotation parameters and displacement parameters must be calculated in turn, but the calculation may cause error accumulation in certain cases. 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. The representative work on dual quaternions was first introduced by Walker et al. (1991). However, the scale parameter in the process of transformation was not considered. Therefore, the scale parameter cannot be directly introduced into the coordinate transformation in geodesy, absolute orientation in photogrammetry, and registration of LiDAR points, which are the main focus of this paper. To compensate for this deficiency, we introduce the scale parameter in the process of similarity transformation and provide the formula derivation process in detail. We also conduct a comparative analysis of the new algorithm based on two experiments. The remainder of the paper is organized as follows. Section 2 briefly reviews the representative work of iterative methods and closed-form methods. Section 3 introduces the concept and properties of dual quaternions and then derives the optimization function based on the Bursa model. The step-by-step and detailed derivation of the formula is also provided. In Section 4, two experiments are described: the first one verifies the correctness of the new algorithm and the second one verifies its feasibility in the registration of LiDAR point clouds. Section 5 concludes the paper.

2. Literature review Most available solutions to the pairwise registration are based on Chasles’ theorem, which states that any rotation and translation can be expressed as a translation along a line and a rotation around

that line (Chasles, 1830). Accordingly, available solution to the pairwise registration problem can be categorized into two, namely, iterative methods and closed-form methods. The iterative method usually uses the rotation matrix to represent rigid transformation and needs the initial approximate values of transformation for linearization, which often introduces error in calculation. By contrast, the closed-form method does not need initial approximate values and avoids linearization, which has been the focus of considerable attention recently. 2.1. Closed-form solution to the pairwise registration problem The most important advantage of the closed-form solution is that it provides the best possible transformation in one step. In contrast to the iterative approach, the closed-form solution does not need to determine good initial estimates. Representatives of closed-form solutions to the pairwise registration problem were independently proposed by Arun et al. (1987), Horn (1987,1988), and Haralick et al. (1989), which differ with respect to transformation representation and alternative ways of minimizing a criterion function. Arun et al. (1987) used the rotation matrix to represent spatial transformation and determined the least squares solution by computing the SVD of the derived matrix (Haralick et al., 1989). Similarly, Horn used orthonormal matrices to represent rotation. After the computation of the positive semidefinite square root of a positive semidefinite matrix, a solution for the rotation matrix is derived through the direct manipulation of 3  3 matrices (Horn, 1987). The third study, which is also similar to that of Horn, used unit quaternions to represent rotations. The solutions for the desired quaternion were the eigenvector of a symmetric 4  4 matrix associated with the largest positive eigenvalue (Horn, 1988). The fourth study is that of Walker et al. (1991), which used dual quaternions to represent transform components. By minimizing a single-cost function associated with the sum of the orientation and position errors, the rotation matrix and the translation vector are derived simultaneously. These closed-form pairwise techniques are theoretically equivalent, but they differ in the way the problem is formulated. An experimental comparison found that the differences were negligible in practical applications with nondegenerate data (Eggert et al., 1997). 2.2. Quaternion-based registration methods Quaternions were invented by W.R. Hamilton to represent a 3D vector (Hamilton, 1844), which is a convenient way to describe rotation on a unit sphere. After Horn successfully applied quaternions in absolute orientation (Horn, 1987), their compactness in describing the rotation matrix and the high efficiency of this description became the focus of considerable attention. Currently, quaternions have been successfully used in rigid body motion analysis (Joseph and Laviola, 2003; Kim and Golnaraghi, 2004) and datum transformation in geodesy (Yang, 1999; Shen et al., 2006; Zeng and Yi, 2011; Li et al., 2012). However, with unit quaternions representing rotation in space, seven transformation parameters are calculated in turn, that is, rotation angles first, scale parameters later, and translation vectors last. In the case of introducing error in calculating the rotation angles, the error may be amplified using the following two steps. To overcome the aforementioned limitation, Walker et al. (1991) introduced dual quaternions to represent similarity transformation. A single-cost function is formulated, with the real part and the dual part of the dual quaternion used to represent rotation and translation, which enable the simultaneous calculation of six parameters (Walker et al., 1991), namely, three rotation angles and three translation parameters. A similar work was conducted by Daniilidis (1999). However, the studies

65

Y. Wang et al. / ISPRS Journal of Photogrammetry and Remote Sensing 94 (2014) 63–69

of Walker et al. (1991) and Daniilidis (1999) do not consider the scale parameter in the process of transformation. 3. Dual quaternion-based registration model

A dual quaternion has eight elements, whereas a 3D object transformation has six independent variables, indicating that two of the eight elements in the dual quaternion representation are not independent. As shown in Eqs. (5) and (6), the components of any dual quaternion satisfy the following two constraints:

3.1. Dual quaternions

r_ T r_ ¼ 1;

ð7Þ

Dual quaternions are composed of quaternions and dual numbers expressed as follows:

r_ T s_ ¼ 0:

ð8Þ

_

q ¼ r_ þ es_ ;

ð1Þ

where r_ and s_ are both real quaternions and are called the real part and the dual part, respectively. Dual quaternions are similarly interpreted as real quaternions:

2

_ 3 cos 2h 6 7 q¼4 _ _ 5; sin 2h n

_

ð2Þ _

where the dual vector n represents a line in 3D space_from which the coordinate system is rotated and translated, and h is the dual angle of rotation and translation. Detailed descriptions of the dual _ _ vector n and dual angle h are expressed as follows: _

n ¼~ n þ e~ p ~ n;

The feature points of the reference station and unregistered station are represented by f~ p0i g and f~ p0i g, respectively, and integer NðN P 3Þ denotes the number of corresponding feature point pairs. Based on the similarities in the feature points between the reference station and the unregistered station, the error norm for the similarity measurement of the proposed registration algorithm can be expressed as follows:

Fðl; R;~ tÞ ¼

ð4Þ

where ~ n is a unit vector that specifies the direction of the rotation axis and the direction of the translation, with the rotation involving the line with direction ~ n passing through point P with a rotation angle of h, and d is the distance of translation along the direction specified by ~ n. Comparing the unit quaternions with the dual quaternion representation, the same transformation can be formed by first translating the original coordinate frame along the direction of ~ n by a distance of d and then rotating it by an angle of h with respect to a line having a unit vector ~ n as its direction and passing through point P. Substituting Eqs. (3) and (4) into Eq. (2), we obtain the following equations (see Fig. 1):

" r_ ¼

# cos 2h h  ; sin 2 ~ n

"

ð5Þ  h

d 2

  1 0 : 2 ~ p

ð10Þ

Based on Eq. (5), the real part of dual quaternions fully represents the whole process of rotation, and the corresponding rotation matrix is expressed as follows:

R ¼ ðr 20  ~ r T~ rÞI þ 2~ r~ rT þ 2r 0 Kð~ rÞ; 2

ð11Þ

3

0 r 3 r2 where Kð~ 0 r1 5, which is a skew-symmetric rÞ ¼ 4 r 3 r2 r 1 0 matrix. Eq. (11) can be expressed as a matrix:

"

R ~ 0T

~ 0 1

# T

¼ Wðr_ Þ Q ðr_ Þ;

where Wðr_ Þ ¼



r0 ~ r

ð12Þ

~ rT r0 I  Kð~ rÞ



and Q ðr_ Þ ¼



r0 ~ r

 ~ rT . ~ r0 I þ KðrÞ

t can be expressed with dual quaternions as Translation vector ~ follows:

#

 d sin 2 : h 2  cos 2 ~ p ~ nÞ n þ sin 2h ð~

ð9Þ

p is The location quaternion of a feature point ~

p_ ¼

h ¼ h þ ed;

N X t ~ klR~ pi þ ~ p0i k: i¼1

ð3Þ

_

s_ ¼

3.2. Error norm and its minimization

ð6Þ

T t_ ¼ Wðr_ Þ s_ ;

ð13Þ  1



t. where t_ T ¼ 2 0 ~ t T , which is the corresponding form of vector ~ Assuming that p_ 0i is the location quaternion of unregistered feature points and p_ i is the new location of quaternion after registration, the dual quaternion-based registration process can be expressed as follows: T T p_ i ¼ Wðr_ Þ s_ þ lWðr_ Þ Q ðr_ Þp_ 0i :

ð14Þ

Therefore, the sum of the square distances between p_ i and p_ 0 is 2

fi ¼ ðp_ i  p_ 0i Þ :

ð15Þ

Substituting Eq. (14) into Eq. (15), we derive the following equation: 2

2 T T fi ¼ ðp_ i  p_ 0i Þ ¼ ðWðr_ Þ s_ þ lWðr_ Þ Q ðr_ Þp_ 0i  p_ 0i Þ 2

2

T

T T T T T ¼ ½Wðr_ Þ s_  þ ½lWðr_ Þ Q ðr_ Þp_ 0i  þ ðp_ 0i Þ ðp_ 0i Þ þ 2l½Wðr_ Þ s_  ½Wðr_ Þ Q ðr_ Þp_ 0i  : T T T T 2ðp_ 0i Þ ½Wðr_ Þ s_   2lðp_ 0i Þ ½Wðr_ Þ Q ðr_ Þp_ 0i 

Fig. 1. Dual quaternion-based rotation and translation.

ð16Þ

66

Y. Wang et al. / ISPRS Journal of Photogrammetry and Remote Sensing 94 (2014) 63–69

Then, each item of Eq. (16) can be further expressed as follows: 2

T T T T ½Wðr_ Þ s_  ¼ ½Wðr_ Þ s_ ½Wðr_ Þ s_  ¼ Wðr_ Þ Wðr_ Þs_ T s_ ¼ ½r_ T r_ s_ T s_

¼ s_ T s_ ;

ð17Þ

r_ T C1 r_ þ s_ T C2 r_ : 2C 4

l¼

T

_0 ¼ l2 ½Wðr_ Þ Wðr_ Þ½Q ðr_ Þ Q ðr_ Þ½p_ 0T i pi  2

_T _

_T _

¼ l ½r r I½r

_0 rI½p_ 0T i pi 

¼l

2 _ 0T _ 0 pi pi ;

ð18Þ

T T T T T 2ðp_ 0i Þ ½Wðr_ Þ s_  ¼ 2ðp_ 0i Þ Wðr_ Þ s_ ¼ 2ðWðr_ Þp_ 0i Þ

¼

T 2ðQ ðp_ 0i Þr_ Þ s_

_T

¼ 2r

T 2ðp_ 0i Þ ½

lWðr_ Þ

T

Q ðr_ Þp_ 0i 

¼ 2l

T ðWðr_ Þp_ 0i Þ ½Wðp_ 0i Þr_ 

¼ 2lr

ð20Þ

2 fi ¼ ðp_ i  p_ 0i Þ

ð32Þ



 1 T T C 22 : r_ C3 C2 r_ = 2C 4  2N 2N

ð22Þ

Using Eqs. (7) and (8) as restrictions, the best dual quaternions for representing rotation can be calculated by minimizing following error function: F ¼ lr_ T C1 r_ þ Ns_ T s_ þ ls_ T C2 r_ þ s_ T C3 r_ þ l2 C 4 þ C 5 þ k1 ðr_ T r_  1Þ þ k2 ðs_ T r_ Þ; ð23Þ

@F ¼ @ r_



lðC1 þ CT1 Þ  lCT2



  

1 1 ðlC2 þ C3 Þ  CT3 ðlC2 þ C3 Þ þ 2k1 r_ 2N 2N ð34Þ

As C3 is skew-symmetric, we derive the following equation:

CT3 C3 ¼ C 33 E:

ð35Þ

Using Eqs. (32) and (35), Eq. (34) can be simplified as follows:

@F 1 2 1 1 1 ¼ lðC1 þ CT1 Þ  l C 22 E  l CT2 C3  l CT3 C2  C 33 E þ 2k1 r_ @ r_ 2N 2N 2N 2N

1 2 1 1 T T C 33 E  lC2 C3 þ lðC1 þ C1 Þ þ 2k1 r_ ¼  l C 22 E  2N 2N N  

1 l 1 T 2 ðl C 22 þ C 33 ÞE  C2 C3  ðC1 þ CT1 Þ þ k1 r_ ¼ 0; ð36Þ ¼2  4N 2 N making

  1 l 1 T ðl2 C 22 þ C 33 ÞE þ C2 C3  ðC1 þ CT1 Þ : 4N 2 N

where k1 and k2 are Lagrange multiplier constants. Taking the partial derivatives, we derive the following equations:

A¼

@F ¼ lðC1 þ CT1 Þr_ þ lCT2 s_ þ CT3 s_ þ 2k1 r_ þ k2 s_ ¼ 0; @ r_

ð24Þ

Eq. (36) can be rewritten as

@F ¼ 2Ns_ þ lC2 r_ þ C3 r_ þ k2 r_ ¼ 0; @ s_

ð25Þ

@F ¼ r_ T C1 r_ þ s_ T C2 r_ þ 2lC 4 ¼ 0: @l

ð26Þ

Therefore, the solution of Eqs. (7), (8), (24), (25), and (26) for r_ and s_ provides the optimal solution for the registration parameters. Multiplying Eq. (25) with r_ T , we obtain the following equation:

k2 ¼ r_ T ðlC2 r_ þ C3 r_ Þ:

ð27Þ

Considering that C2 and C3 are both skew-symmetric matrices, we obtain the following equation:

k2 ¼ r_ T ðlC2 r_ þ C3 r_ Þ ¼ 0:

ð33Þ

Substituting Eq. (29) into Eq. (24), we obtain the following equation:

ð21Þ

P P P T making C1 ¼ 2 Ni¼1 Q ðp_ 0i Þ Wðp_ 0i Þ, C2 ¼ 2 Ni¼1 Wðp_ 0i Þ, C3 ¼ 2 Ni¼1 P P ~0 Q ðp_ 0i Þ, C 4 ¼ Ni¼1 ð~ pTi ~ pi Þ, and C 5 ¼ Ni¼1 ð~ p0T i pi Þ. The similarity measurement, as shown in Eq. (9), can be rewritten in the form of the quadratic function of quaternions r_ and s_ :

F ¼ lr_ T C1 r_ þ Ns_ T s_ þ ls_ T C2 r_ þ s_ T C3 r_ þ l2 C 4 þ C 5 :

ð31Þ

¼ 0:

T

_0 _0 _0 _T _0 _ _T _ 0 _ ¼ s_ T s_ þ l2 p_ 0T i pi þ ðpi Þ ðpi Þ þ 2ls Wðpi Þr  2s Q ðpi Þr T

:

where E is a unit matrix. Using Eqs. (32), Eq. (31) can be further expressed as



Substituting Eqs. (17)–(20) into Eq. (16), we derive the following equation:

 2lr_ T Q ðp_ 0i Þ Wðp_ 0i Þr_ ;

2C 4

l ¼ r_ T C1 r_ þ

T

¼ 2lðQ ðp_ 0i Þr_ Þ ½Wðp_ 0i Þr_  T Q ðp_ 0i Þ Wðp_ 0i Þr_ :

l¼

CT2 C2 ¼ C 22 E; ð19Þ

_T

h i l _T T 1 _T T r_ T C1 r_  2N r C2 C2 r_ þ 2N r C3 C2 r_

As C2 is skew-symmetric, we derive the following equation:

T Q ðp_ 0i Þ s_

¼ 2s_ T Q ðp_ 0i Þr_ ;

ð30Þ

Substituting Eq. (29) into Eq. (30), we obtain the following equation:

2

T T T ½lWðr_ Þ Q ðr_ Þp_ 0i  ¼ l2 ½Wðr_ Þ Q ðr_ Þp_ 0i ½Wðr_ Þ Q ðr_ Þp_ 0i  T

Using Eq. (26), scale parameter l can be represented as the function of r_ and s_ as follows:

ð28Þ

Ar_ ¼ k1 r_ :

ð37Þ

ð38Þ

Based on Eq. (38) and according to Walker et al. (1991), the quaternion r_ is an eigenvector of the matrix A and k1 is the corresponding eigenvalue. In general, four solutions are derived for this equation. Considering that A is real and symmetric, all the eigenvalues and eigenvectors are real and all the eigenvectors are orthogonal. The desired solution is identified by referring back to the original error shown in Eq. (22). Multiplying Eq. (24) by r_ T , we obtain the following equation:

1 1 lr_ T C1 r_ ¼  lr_ T CT2 s_  r_ T CT3 s_  k1 : 2

2

ð39Þ

Then, multiplying Eq. (25) by s_ T , we obtain:

1 Ns_ T s_ ¼  ðls_ T C2 r_ þ s_ T C3 r_ Þ: 2

ð40Þ

Substituting Eqs. (39) and (40) into Eq. (22) gives

Therefore, we can represent the dual quaternion s_ as the function of r_ as follows:

F ¼ l2 C 4 þ C 5  k1 :

1 ðlC2 r_ þ C3 r_ Þ: s_ ¼  2N

Clearly, the error will be minimized if we select the eigenvector corresponding to the largest positive eigenvalue.

ð29Þ

ð41Þ

67

Y. Wang et al. / ISPRS Journal of Photogrammetry and Remote Sensing 94 (2014) 63–69 Table 1 Coordinate pairs extracted from the local and WGS84 coordinate systems. Site

Solitude Buoch Zeil Hohenneuffen Kuehlenberg Ex Mergelaec Ex Hof Asperg Ex Kaisersbach

Local coordinate system (System A)

WGS84 coordinate system (System B)

x

y

z

X

Y

Z

4157222.543 4149043.336 4172803.511 4177148.376 4137012.190 4146292.729 4138759.902

664789.307 688836.443 690340.078 642997.635 671808.029 666952.887 702670.738

4774952.099 4778632.188 4758129.701 4760764.800 4791128.215 4783859.856 4785552.196

4157870.237 4149691.049 4173451.354 4177796.064 4137659.549 4146940.228 4139407.506

664818.678 688865.785 690369.375 643026.700 671837.337 666982.151 702700.227

4775416.524 4779096.588 4758594.075 4761228.899 4791592.531 4784324.099 4786016.645

Table 2 Contrast of the proposed dual quaternion-based algorithm and the weighted Procrustes algorithm. Scheme No. Dual quaternion-based algorithm

Weighted Procrustes algorithm

Rotation matrix (R) 2 3 1:000000000 0:000004815 0:000004333 4 0:000004815 1:000000000 0:000004841 5 0:000004333 0:000004841 1:000000000 2 3 0:999999999 0:000004814 0:000004332 4 0:000004814 0:999999999 0:000004841 5 0:000004332 0:000004841 0:999999999

Translation vector (T)

Scale (l)

(641.8805, 68.6554, 416.3982)

1.000005583

(641.8804, 68.6553, 416.3982)

1.000005583

Table 3 Residuals of each site after transformation. Site

Solitude Buoch Zeil Hohenneuffen Kuehlenberg Ex Mergelaec Ex Hof Asperg Ex Kaisersbach

Dual quaternion-based algorithm

Weighted Procrustes algorithm

Dx

Dy

Dz

DX

DY

DZ

0.0939 0.0588 0.0403 0.0197 0.0916 0.0117 0.0292

0.1353 0.0500 0.0883 0.0213 0.0140 0.0067 0.0035

0.1402 0.0136 0.0078 0.0872 0.0059 0.0549 0.0014

0.0940 0.0588 0.0399 0.0202 0.0919 0.0118 0.0294

0.1351 0.0497 0.0879 0.0220 0.0139 0.0065 0.0041

0.1402 0.0137 0.0081 0.0874 0.0055 0.0546 0.0017

  Making A0 ¼ N1 CT2 C3  C1 þ CT1 and kA0 as its eigenvector, the eigenvalue of matrix A can be expressed as

1 l kA ¼ ðl2 C 22 þ C 33 Þ þ kA0 : 4N 2

Table 4 Point features extracted from two neighboring LiDAR point clouds. No.

ð42Þ

Considering that l, N, C22, and C33 are all constants, kA is maximized only when kA00 is maximized. Therefore, the calculation of the maximum eigenvalue of matrix A0 and its eigenvector r_ efficiently obtains our desired result By substituting r_ into Eq. (31), (33), l is obtained. Using Eq. (29) results in s_ . Therefore, the best registration parameters are acquired. 3.3. Implementation of the Proposed Registration Algorithm The detailed implementation of the proposed registration algorithm can be summarized as follows: 1. We input the two sets of corresponding feature points f~ p0i g and 0 f~ pi g, which are extracted from the reference station and unregistered station, respectively. 2. We construct matrix A0 based on matrices C1, C2, and C3 as well as the two constants C4 and C5. Then, we calculate the maximum eigenvalue and its corresponding eigenvector r_ of A0 . 3. We calculate the scale parameter l using Formula (31), or (33) and the quaternion s_ using Formula (29). 4. We calculate the rotation matrix R and translation vector ~ t using Formula (13).

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

Reference station

Unregistered station

x

y

z

x

y

z

91.406 91.297 60.158 60.135 56.298 13.269 4.666 49.939 52.769 72.929 46.500 52.581 58.972 55.429 55.313 63.467 57.673 49.687

53.344 53.222 24.280 24.278 19.186 2.677 17.245 14.297 11.523 8.630 30.291 22.934 17.511 26.155 26.131 27.962 22.069 14.083

8.320 0.916 8.948 1.521 5.700 1.444 1.605 27.119 25.906 27.146 23.078 5.676 18.862 23.077 23.039 26.981 25.782 3.666

49.007 47.365 36.514 34.881 53.378 7.324 9.587 36.532 39.932 67.051 54.124 51.943 57.712 59.650 59.512 41.466 39.133 29.781

54.453 54.435 13.733 13.859 25.872 32.695 19.650 0.319 1.307 8.834 40.688 30.962 23.376 32.625 32.705 18.246 10.234 0.026

0.978 6.242 3.642 3.608 4.187 1.389 2.449 21.980 19.965 15.017 13.216 3.965 8.397 12.037 12.071 21.085 20.247 8.062

4. Experiments 4.1. Experimental scheme designs The proposed algorithm is implemented using C++ programming language. Two experimental schemes are designed to verify the correctness and effectiveness of the algorithm. The first

68

Y. Wang et al. / ISPRS Journal of Photogrammetry and Remote Sensing 94 (2014) 63–69

Table 5 Results of the proposed dual quaternion-based registration algorithm. Scheme No. 1

2

Rotation matrix (R) 2 0:850416493 0:494507078 4 0:479380907 0:868981201 0:216761931 0:018287246 2 0:850416493 0:494507078 4 0:479380907 0:868981201 0:216761931 0:018287246

3

0:179595486 0:122742085 5 0:976053196 3 0:179595486 0:122742085 5 0:976053196

experiment involves the coordinate transformation between two different coordinate systems. The results are compared with those of the well-known weighted Procrustes algorithm (Grafarend and Awange, 2003). The second experiment involves the use of the proposed algorithm in the registration of LiDAR points. The correctness of the proposed model after the introduction of the scale parameter is analyzed in detail. 4.2. Coordinate transformation in geodesy The Cartesian coordinates of the seven stations, as shown in Table 1, are taken from Grafarend and Awange (2003). Using these coordinates, we compute the transformation parameters from a local geodetic system (System A) to WGS84 (System B). The transformation parameters from the local coordinate system to WGS84 are calculated by the proposed dual quaternion-based algorithm and the weighted Procrustes algorithm. The results are used to verify the correctness of the dual quaternion-based algorithm. The results of the two algorithms are shown in Table 2. The residuals of each site after transformation are shown in Table 3. Based on the residuals shown in Table 3, the standard error of the dual quaternion-based algorithm is 0.011799 and that of the weighted Procrustes algorithm is 0.011796. The difference is at the fifth place after the decimal point, which is small enough to ignore. Therefore, we can conclude that the dual quaternion-based algorithm is valid to calculate the similarity transformation parameters. 4.3. LiDAR point cloud registration The second experimental data are the two LiDAR point clouds acquired by a terrestrial LiDAR instrument (Riegl LMS-Z420i). Currently, most registration algorithms are highly efficient, and almost all of them can calculate real-time or quasi-real-time transformation parameters. Therefore, in this study, only the correctness of the presented algorithm is verified. Sample point clouds are acquired by the Riegl LMS-Z420i series terrestrial LiDAR. Conjugate point features are extracted by (1) fitting salient points marked by special reflective material and (2) by taking the intersection of three planes. Both extraction means are highly accurate. The extracted conjugate point features used to test the proposed algorithm are shown in Table 4. The extracted conjugate point features are directly included in the presented dual quaternion-based registration algorithm. The transformation parameters are presented in the first row in Table 5. To further verify the validity of the presented algorithm, all point coordinates from the unregistered station are scaled to half of their actual coordinates and are included in the proposed dual quaternion-based registration algorithm. The transformation parameters are shown in the second row in Table 5. To further verify the correctness of the dual quaternion-based registration algorithm, the unit quaternion-based algorithm (Horn, 1987) is used to calculate the transformation parameters. Comparatively, the results of the unit quaternion-based algorithm are the same as those in Table 5.

Translation vector (T)

Scale factor (s)

22.9656, 29.3962, 2.2652

1.000385435

22.9656, 29.3962, 2.2652

2.000770870

Therefore, we conclude that the proposed dual quaternionbased registration algorithm works and is an alternative for the concise description of similarity transformation. 5. Conclusions The representation of rigid transformation in 3D space, especially of orientation, is a crucial issue in navigation, geodesy, photogrammetry, and robot arm manipulation, among others. Currently, the most popular representations of a 3D rotation are rotation matrix, Euler angles, Rodrigues vector, and quaternions. As renormalization is difficult for rotation matrices, representation by Euler angles may cause singularities, and Rodrigues vectors do not allow for an easy composition algorithm, quaternions seem to be a suitable representation of rotations in 3D with few parameters. Considering that unit quaternions represent rotation in space, seven transformation parameters are calculated, namely, the rotation angles first, the scale parameter later, and three translation vectors last. In the case of introducing an error in calculating the rotation angles, the error may be amplified in the following two steps. Dual quaternions are used to represent similarity transformation. A single-cost function is formulated in which the real part and the dual part of the dual quaternions are used to represent rotation and translation. By minimizing the formulated single-cost function, seven transformation parameters are calculated simultaneously. Detailed derivations of related formulas are given and analyzed, and experiments are conducted to verify the effectiveness and the feasibility of the proposed algorithm. As only point features are used to construct the similarities between the two neighboring LiDAR stations, future work can focus on introducing line and face features in pairwise registration. A newly formulated dual quaternion-based error norm should consider orientation and position errors simultaneously to formulate a unified pairwise registration model that considers point, line, and face features simultaneously. Acknowledgments This study was funded by the National Natural Science Foundation of China (Project Numbers 41271444/41001297/41001312) and the Natural Science Foundation of Jiangsu Province (Project Number BK2012569). References Arun, K.S., Huang, T.S., Blostein, S.D., 1987. Least-squares fitting of two 3-D point sets. IEEE Trans. Pattern Anal. Mach. Intell. 9 (5), 698–700. Besl, P.J., McKay, N.D., 1992. A method for registration of 3-D shapes. IEEE Trans. Pattern Anal. Mach. Intell. 14 (2), 239–256. Chasles, M., 1830. Note sur les propriétés générales du système de deux corps semblables entr’eux et placés d’une manière quelconque dans l’espace; et sur le déplacement fini ou infiniment petit d’un corps solide libre. Bull. Sci. Math. Astron. Phys. Chim. 14 (XIV), 321–326. Daniilidis, K., 1999. Hand-eye calibration using dual quaternions. Int. J. Robotics Res. 18 (3), 286–298. Eggert, D.W., Lorusso, A., Fisher, R.B., 1997. Estimating 3-D rigid body transformations: a comparison of four major algorithms. Mach. Vis. Appl. 9 (5–6), 272–290.

Y. Wang et al. / ISPRS Journal of Photogrammetry and Remote Sensing 94 (2014) 63–69 Grafarend, E.W., Awange, L.J., 2003. Nonlinear analysis of the three dimensional datum transformation [conformal group C7(3)]. J. Geodesy 77 (1), 66–76. Hamilton, W.R., 1844. On quaternions, or on a new system of imaginaries in algebra. Phil. Mag. 25 (3), 489–495. Haralick, R.M., Joo, H., Lee, C.N., et al., 1989. Pose estimation from corresponding point data. IEEE Trans. SMC 19 (6), 1426–1446. Horn, B.K., 1987. Closed-form solution of absolute orientation using unit quaternions. J. Opt. Soc. Am. (Ser. A) 4 (4), 629–642. Horn, B.K., 1988. Closed form solution of absolute orientation using orthonormal matrices. J. Opt. Soc. Am. (Ser. A) 5 (7), 1127–1135. Prošková, J., 2011. Application of dual quaternions algorithm for geodetic datum transformation. J. Appl. Math. 4 (2), 225–236. Prošková, J., 2012. Discovery of dual quaternions for geodesy. J. Geometry Graphics 16 (2), 195–209. Joseph, J., Laviola, J., 2003. A comparison of unscented and extended Kalman filtering for estimating quaternions motion. IEEE Press. 3, 2435–2440.

69

Kim, A., Golnaraghi, M.F., 2004. A quaternions-based orientation estimation algorithm using an inertial measurement unit. IEEE Position Location Navigation Symp., 268–272. Li, B., Shen, Y., Li, W., 2012. The seamless model for three-dimensional datum transformation. Sci. China (Earth Sci.) 55 (12), 2099–2108. Walker, M.W., Shao, L., Volz, R.A., 1991. Estimating 3-D location parameters using dual number quaternions. CVGIP: Image Understanding 54 (3), 358–367. Yang, Y., 1999. Robust estimation of geodetic datum transformation. J. Geodesy 73 (5), 268–274. Shen, Y.Z., Chen, Y., Zheng, D.H., 2006. A quaternions-based geodetic datum transformation algorithm. J. Geodesy 80 (5), 233–239. Zeng, H., Yi, Q., 2011. Quaternion-based iterative solution of three-dimensional coordinate transformation problem. J. Comput. 6 (7), 1361–1368. Zheng, D., Yue, D., Yue, J., 2008. Geometric feature constraint based algorithm for building scanning point cloud registration. Acta Geodaetica Et Cartographica Sin. 37 (4), 464–468.