1 Introduction
We consider the problem of single rotation averaging, i.e., averaging several estimates of a single rotation to obtain the best estimate. This problem is relevant in many applications such as structure from motion (SfM)
[7, 13], camera rig calibration [3], motion capture [12], satellite/spacecraft attitude determination [9, 10] and crystallography [8, 11].A standard approach for single rotation averaging is to find the rotation that minimizes a cost function based on the distance to the input rotations. We refer to [6] for an extensive study of various distance functions. The current stateoftheart method is to minimize the sum of geodesic distances using the Weiszfeld algorithm on [7].
In this work, we propose a novel method, also based on the Weiszfeld algorithm [14, 15], that is faster and more robust than [7]. Our contributions are as follows:

A robust initialization from the elementwise median of the input rotation matrices (Section 3.1).

An implicit outlier rejection scheme performed at each iteration of the Weiszfeld algorithm (Section 3.2).

An approximation of the chordal median in using the Weiszfeld algorithm (Section 3.3).
We substantiate our claim through extensive evaluation on synthetic data (Section 4).
2 Preliminaries
We denote the vectorization of an
matrix by and its inverse by . For a 3D vector , we define as the corresponding skewsymmetric matrix, and denote the inverse operator by , i.e., . The Euclidean, the and the Frobenius norm are respectively denoted by , and . A rotation can be presented by a rotation matrix or a rotation vector where and represent the angle and the unit axis of the rotation, respectively. The two representations are related by Rodrigues formula, and we denote the corresponding mapping between them by and [5]:(1) 
(2) 
(3) 
The geodesic distance between two rotations is obtained by substituting into (3). In [6], the relation between the geodesic and chordal distance was given:
(4)  
(5) 
We define as the projection of the matrix onto the special orthogonal group , which gives the closest rotation in the Frobenius norm [2]: For ,
(6) 
where
(7)  
(8) 
3 Method
3.1 Robust Initialization
In [7], the chordal mean of the rotations is taken as the starting point of the Weiszfeld algorithm. For input rotations , it is given by [6]. Although this initial solution can be obtained very fast, it is often inaccurate and sensitive to outliers. To overcome this weakness, we propose to initialize using the following matrix:
(9) 
where the subscript denote the element at the th row and the th column of the matrix. Note that is called the elementwise norm of the matrix . See Fig. 1 for the geometric interpretation of this distance metric. Since the nine entries of are independent, we can consider them separately in 1D space. Then, the entry of at location minimizes the sum of absolute deviations from the entries of ’s at , meaning that it is simply their median:
(10) 
The initial rotation matrix is then be obtained by projecting onto :
(11) 
3.2 Outlier Rejection in the Weiszfeld Algorithm
The geodesic mean (i.e., median) of the rotations is defined as
(12) 
In [7], it was shown that this can be computed using the Weiszfeld algorithm on and that it is more robust to outliers than the mean. However, a large number of outliers is still critical to the accuracy. To further mitigate the influence of the outliers, we modify the Weiszfeld algorithm of [7] such that the large residuals are given zero weight at each iteration. Specifically, we disregard all the residuals larger than where
is the first quartile of the residuals at each iteration and
is some threshold we set in order to avoid discarding inliers. The details are given in Algorithm 1. A similar approach of disregarding large residuals was used in [4] for robust PerspectivePoint (PnP) problem. Note that our approach contrasts [1] where a smooth robust cost function is favored for theoretically guaranteed convergence. In practice, our method is more robust to outliers (see Section 4.2).3.3 Approximate Chordal Mean
The chordal mean of the rotations is defined as
(13) 
In [6], a locally convergent algorithm on is proposed for this problem. In this work, we propose a different approach: Instead of iteratively updating the estimate on , we first embed the rotations in a Euclidean space , find their geometric median in using the standard Weiszfeld algorithm [14, 15] (which is globally convergent), and then project this median onto . In other words, we approximate as
(14) 
(15)  
(16) 
Since the optimization is performed using the Weiszfeld algorithm, we can also incorporate the initialization and outlier rejection scheme in the previous sections. Algorithm 2 summarizes our method.
We point out two things in the implementation: First, since we do not optimize on , the initial estimate does not have to come from a rotation, and we omit (11). Second, the threshold must be scaled appropriately when comparing Algorithm 1 and 2. Assuming that at each iteration does not vastly differ from an embedding of a rotation in , we convert from geodesic to chordal using (5), and vice versa. This is done in line 2 of Algorithm 2.
4 Results
4.1 Initialization
For evaluation, we generated a synthetic dataset where the inlier rotations follow a Gaussian distribution with
, and the outliers have uniformly distributed angles
at random directions. Fig. 2 compares the average accuracy of the proposed initial solution (Section 3.1) and the chordal mean [6] over 1000 runs. It can be seen that our solution is significantly better than the chordal mean unless the outlier ratio is extremely high (i.e., above ). On average, the chordal method takes 0.37 s and ours 0.83 s per rotation. This time difference is insignificant compared to the optimization that follows (see Tab. 1).4.2 Comparison against [7]
Using the same setup as in previous section, we compare Algorithm 1 and 2, with and without the proposed outlier rejection scheme (Section 3.2). This time, we consider two different inlier noise levels, and . The average accuracy of the evaluated methods^{1}^{1}1We did not include the chordal mean here, since it produced much larger errors than the rest and was already reported in Fig. 2. is compared in Fig. 3. With the outlier rejection, the geodesic mean and our approximate chordal mean are almost equally accurate. Without the outlier rejection, the geodesic mean is more accurate than our approximate chordal mean, but only for very high outlier ratios (i.e, ). Otherwise, there is no significant difference between the two.
The computation times are reported in Tab. 1. Our method is always faster than [7], and is 2–4 times faster with the outlier rejection. That said, the speed is not a major advantage, since all methods can process several hundreds of rotations in less than a millisecond. In most cases, averaging rotations will take much less time than other operations, such as the computation of input rotations.
5 Conclusions
In this work, we proposed a novel alternative to the work of Hartley et al. [7] for robust single rotation averaging. While both our method and [7] use the Weiszfeld algorithm, there are three key differences:

We initialize the Weiszfeld algorithm using the elementwise median of the input rotation matrices.

We implicitly disregard the outliers at each iteration of the Weiszfeld algorithm.

We approximate the chordal median on instead of the geodesic median as in [7].
As a result, our method achieves better performance in terms of speed and robustness to outliers. We also found that incorporating the proposed outlier rejection in the original implementation of [7] leads to similar performance, but at 2–4 times slower speed than ours.
w/o outlier rejection  w/ outlier rejection  

[7]  Ours  [7]  Ours  
(, )  8.69  4.38 (2.0)  9.37  4.43 (2.1) 
(, )  10.5  4.47 (2.3)  11.7  5.68 (2.1) 
(, )  15.2  6.21 (2.4)  17.2  6.78 (2.5) 
(, )  24.9  15.2 (1.6)  27.2  7.71 (3.5) 
(, )  32.1  10.3 (3.1)  31.6  8.99 (3.5) 
(, )  10.7  6.00 (1.8)  17.0  6.02 (2.8) 
(, )  15.0  5.98 (2.5)  17.0  7.1 (2.4) 
(, )  19.6  7.66 (2.6)  22.3  8.06 (2.8) 
(, )  24.9  11.1 (2.2)  28.1  8.77 (3.2) 
(, )  29.1  10.2 (2.9)  31.7  8.50 (3.7) 
References

[1]
K. Aftab and R. Hartley.
Convergence of iteratively reweighted least squares to robust
mestimators.
In
IEEE Winter Conf. on Applications of Computer Vision
, pages 480–487, 2015.  [2] K. S. Arun, T. S. Huang, and S. D. Blostein. Leastsquares fitting of two 3D point sets. IEEE Trans. Pattern Anal. Mach. Intell., 9(5):698–700, 1987.
 [3] Yuchao Dai, Jochen Trumpf, Hongdong Li, Nick Barnes, and Richard Hartley. Rotation averaging with application to camerarig calibration. In Asian Conf. on Computer Vision, pages 335–346, 2009.

[4]
Luis Ferraz, Xavier Binefa, and Francesc MorenoNoguer.
Very fast solution to the pnp problem with algebraic outlier
rejection.
In
Proc. IEEE Computer Society Conf. on Computer Vision and Pattern Recognition (CVPR)
, pages 501–508, 2014.  [5] Christian Forster, Luca Carlone, Frank Dellaert, and Davide Scaramuzza. Onmanifold preintegration for realtime visual–inertial odometry. IEEE Trans. Robot., 33(1):1–21, 2017.
 [6] Richard Hartley, Jochen Trumpf, Yuchao Dai, and Hongdong Li. Rotation averaging. International Journal of Computer Vision, 2013.
 [7] Richard I. Hartley, Khurrum Aftab, and Jochen Trumpf. L1 rotation averaging using the Weiszfeld algorithm. In IEEE Conf. on Computer Vision and Pattern Recognition (CVPR), pages 3041–3048, 2011.
 [8] M. Humbert, N. Gey, J. Muller, and C. Esling. Determination of a Mean Orientation from a Cloud of Orientations. Application to Electron BackScattering Pattern Measurements. Journal of Applied Crystallography, 29(6):662–666, 1996.
 [9] Q. M. Lam and J. L. Crassidis. Precision attitude determination using a multiple model adaptive estimation scheme. In IEEE Aerospace Conference, pages 1–20, 2007.
 [10] Landis Markley, Yang Cheng, John Crassidis, and Yaakov Oshman. Averaging quaternions. Journal of Guidance, Control, and Dynamics, 30:1193–1196, 07 2007.
 [11] A. Morawiec. A note on mean orientation. Journal of Applied Crystallography, 31(5):818–819, 1998.
 [12] Inna Sharf, Alon Wolf, and M.B. Rubin. Arithmetic and geometric solutions for average rigidbody rotation. Mechanism and Machine Theory, 45(9):1239 – 1251, 2010.
 [13] Roberto Tron, Xiaowei Zhou, and Kostas Daniilidis. A survey on rotation optimization in structure from motion. In IEEE Conf. on Computer Vision and Pattern Recognition Workshops, pages 1032–1040, 2016.
 [14] Endre Weiszfeld. Sur le point pour lequel la somme des distances de n points donnés est minimum. Tohoku Mathematical Journal, 43:355–386, 1937.
 [15] Endre Weiszfeld and Frank Plastria. On the point for which the sum of the distances to n given points is minimum. Annals of Operations Research, 167(1):7–41, 2009.
Comments
intrinsicD ∙
Hello,
I have a question regarding Algorithm 1.
In line 6: what does the notation mean? Do you measure whether R_gm is in the set {R_i}? How do you do that?
In line 7: how are the perturbation matrices R_perturb constructed?
An explanation would be much appreciated!
Thank you and Best Regards.
intrinsicD ∙
turns out i figured it out on my own.