Metric Space for Direction Cosine Matrix?

Question

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

Answer 1

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.

Answer 2

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.