Hello,
The following code is an implementation (and some explanations) of the matrix rotation mean as described in the publication:
R. Hartley, J. Trumpf, Y. Dai, and H. Li, “Rotation Averaging” International Journal of Computer Vision 103, pp. 267–305, DOI 10.1007/s11263-012-0601-0, July 2013 http://users.cecs.anu.edu.au/~hongdong/rotationaveraging.pdf
The first function to be defined is the Skew-Symetrixc matrix (as in https://en.wikipedia.org/wiki/Skew-symmetric_matrix)
function cross_ten_x = cross_ten( x )
% Skew symetric matrix permiting to express the cross product: x * .
X = x(1);
Y = x(2);
Z = x(3);
cross_ten_x = [ 0 , -Z , Y ;...
Z , 0 , -X ;...
-Y , X , 0 ];
end
In math, we use the notation: 
This function have a nice property:
This function have a nice property: for any vector v, 
is the rotation matrix around the axis defined by v, and of angle 
.
is the rotation matrix around the axis defined by v, and of angle .Then, we could define the inverse transformation
function v = lie_log( R )
% Lie-Logarithm
% Function such that:
% expm(cross_ten( lie_log(R) )) = R
[V,D] = eig(R);
[theta,i] = max(angle(diag(D)));
[ ~ ,k] = min(abs(diag(D)-1));
v1 = theta*V(:,k);
v2 =-theta*V(:,k);
R1 = expm(cross_ten( v1 ));
R2 = expm(cross_ten( v2 ));
if max(angle(eig(R * R2'))) > max(angle(eig(R * R1')))
v = v1;
else
v = v2;
end
end
Finally, you could process the rotation matrix mean.
Let Rotations, be the cell array contenning the rotations, then, the followinc code is the implementation of the reference you cited .
%% Example of Input
Rotations = { [ 0.5415 -0.6881 0.4831;... R1
0.8392 0.4769 -0.2613;...
-0.0506 0.5469 0.8357];...
[ 0.5422 -0.6857 0.4857;... R2
0.8386 0.4775 -0.2620;...
-0.0522 0.5494 0.8339];...
[ 0.5430 -0.6857 0.4848;... R3
0.8381 0.4784 -0.2620;...
-0.0523 0.5486 0.8344];...
[ 0.5400 -0.6889 0.4835;... R4
0.8401 0.4764 -0.2593;...
-0.0517 0.5462 0.8360];...
[ 0.5429 -0.6861 0.4843;... R5
0.8381 0.4794 -0.2603;...
-0.0535 0.5472 0.8353] }
%% Algorithm for Matrix averaging
nb_it_max = 20;
tol_r = 1e-10; % [1]
number_rotations = length(Rotations);
R = Rotations{1} % First approx of R [1]
for nb_it = 1:nb_it_max % [2]
list_r = nan(3,number_rotations); % [3]
for i = 1:number_rotations
list_r(:,i) = lie_log( R' * Rotations{i} );
end
r = mean(list_r,2);
if norm(r) < tol_r % [4]
break
end
R = R * expm(cross_ten( r )); % Update [7]
end % [8]
if nb_it == nb_it_max
error('the maximum number of iteration where reached')
end
R
The [*] indicate the number of the line, refering to page 18 of http://users.cecs.anu.edu.au/~hongdong/rotationaveraging.pdf
