`Orientation.dot_outer()` scambles axis ordering and can give the wrong answer

[argerlt]

While writing tests for uncovered lines, I noticed the logic in Orientation.dot_outer scrambles axis ordering for 2D or greater arrays or orientations:

import orix.quaternion as oqu
>>> g1 = oqu.Orientation.random(shape=(2,3,4))
>>> g2 =  oqu.Orientation.random(shape=(5,6))
>>> g12 = g1.dot_outer(g2)
>>> print(g12.shape)
(3,4,5,6,2)

Looking deeper, this logic for finding the smallest angle also isn't guaranteed to work for roto-inversion symmetries:

        symmetry = _get_unique_symmetry_elements(self.symmetry, other.symmetry)
        M = other * ~self
        all_dot_products = Rotation(M).dot_outer(symmetry)
        highest_dot_products = np.max(all_dot_products, axis=-1)
        return highest_dot_products

I fixed both issues in #669, but I'm opening a second issue here just to discuss in case someone disagrees with my method.