# Metric Space for Direction Cosine Matrix?

I work in the medical field. Sometime we receives MR images that have been acquired along the same direction, however when looking up the direction cosine matrix, the values are slightly different (up to a certain precision). The scanner only stores the first 6 values of the matrix, the last 3 values are computed using a cross product from the first 2 vectors. Those two direction cosines define the first row and the first column with respect to the patient.

What is the correct metric space (equivalent to euclidean distance?) to compare two different direction cosine matrix ?

Typical example:

import numpy as np
a1,b1=np.array([0.997704,0.0677201,6.10347E-5]),np.array([-0.0673549,0.992439,-0.102604])
a2,b2=np.array([0.997704,0.0677201,6.10347E-5]),np.array([-0.067355,0.992439,-0.102604])
a3,b3=np.array([0.997704,0.0677202,6.10347E-5]),np.array([-0.0673549,0.992439,-0.102604])


With a lot of hand waving, comparing the dot product of the normal give a sense of how close direction cosine matrix are:

>>> c1 = np.cross(a1,b1)
>>> c2 = np.cross(a2,b2)
>>> c3 = np.cross(a3,b3)
>>> np.dot(c1,c2)
0.9999987261098883
>>> np.dot(c2,c3)
0.9999987328817403

• The function $d(x, y) = \arccos(x\cdot y)$ is a metric on the unit sphere, i.e. it's symmetric, positive-definite, and satisfies the triangle inequality. It's basically just the great circle distance between the two points. Is that what you mean by metric? Sep 23 '21 at 16:03
• @DanielShapero, I think that your comment answer the question. Would you post it as such? Sep 24 '21 at 11:29
• @DanielShapero I've updated my question, since it was not very clear I was talking about direction cosine matrix (not just direction cosines of a single vector). Sep 27 '21 at 6:27

I ended up implementing the suggestion by @DanielShapero (or at least how I understood it):

import numpy as np
def distance(a, b): # a & b contains each the first 6 values of direction cosine matrix
ax, ay = a # ax is the first vector, ay is the second vector
bx, by = b # bx is the first vector, by is the second vector
az = np.cross(ax, ay) # ax × ay
bz = np.cross(bx, by) # bx × by
return np.arccos(np.dot(az, bz)) # d(x,y)=arccos(x⋅y)


Using my original example:

a1, b1 = np.array([0.997704, 0.0677201, 6.10347E-5]), np.array([-0.0673549, 0.992439, -0.102604])
a2, b2 = np.array([0.997704, 0.0677201, 6.10347E-5]), np.array([-0.067355, 0.992439, -0.102604])
a3, b3 = np.array([0.997704, 0.0677202, 6.10347E-5]), np.array([-0.0673549, 0.992439, -0.102604])


as well as:

a4, b4 = np.array([1, 0, 0]), np.array([0, 1, 0])
a5, b5 = np.array([0, 1, 0]), np.array([0, 0, 1])
a6, b6 = np.array([0, 0, 1]), np.array([1, 0, 0])


I can verify this give the expected result:

print(distance((a1, b1), (a2, b2)))
print(distance((a2, b2), (a3, b3)))
print(distance((a3, b3), (a1, b1)))

0.0015961769213946786
0.00159192872155101
0.0015961542345642986


Same for:

print(distance((a4, b4), (a5, b5)))
print(distance((a5, b5), (a6, b6)))
print(distance((a6, b6), (a4, b4)))

1.5707963267948966
1.5707963267948966
1.5707963267948966


Technically the expression based on arctan is the only one that is well-conditioned for all angles. So I ended up not using arccos directly.

If you have a reference direction I would suggest using the norm of the difference of each vector with the reference vector. Since the vectors are already normalized this works as a relative norm. Another option is to compute the angle between your direction and the reference one.

• I do not have a reference direction...unless there is a way to compute the "mean" of multiple direction cosines... Sep 23 '21 at 18:26
• You could compute the average of all your direction vectors. Sep 23 '21 at 20:29
• my head hurts when I read "average" for direction vectors, could you briefly point me to an algorithm ? Sep 24 '21 at 6:06