3
$\begingroup$

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
$\endgroup$
3
  • 2
    $\begingroup$ 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? $\endgroup$ Sep 23 '21 at 16:03
  • $\begingroup$ @DanielShapero, I think that your comment answer the question. Would you post it as such? $\endgroup$
    – nicoguaro
    Sep 24 '21 at 11:29
  • $\begingroup$ @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). $\endgroup$
    – malat
    Sep 27 '21 at 6:27
1
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.