The indirection of anonymous functions costs time. So either use the built-in functions with the same names cross, norm and dot, or hard code the functions directly.
Instead of the expensive trick to determine the positions of instabilities in the ASIN and ACOS methods, use a stable method directly:
atan2(norm(cross(N1 x N2)), dot(N1, N2))
Where N1 and N2 are the normalized input vectors.
N1 = bsxfun(@rdivide, a, sqrt(sum(a .* a ,1)))
N2 = bsxfun(@rdivide, b, sqrt(sum(b .* b ,1)))
N1dotN2 = N1(:, 1) .* N2(:, 1) + N1(:, 2) .* N2(:, 2) + N1(:, 3) .* N2(:, 3);
N1xN2 = [(N1(:, 2) .* N2(:, 3) - N1(:, 3) .* N2(:, 2)), ...
(N1(:, 3) .* N2(:, 1) - N1(:, 1) .* N2(:, 3)), ...
(N1(:, 1) .* N2(:, 2) - N1(:, 2) .* N2(:, 1))];
Angle = atan2(sqrt(sum(N1xN2 .* N1xN2, 1)), N1dotN2);
Best Answer