0

I've read through this question (Understanding DICOM image attributes to get axial/coronal/sagittal cuts) and responses a bunch, but I still can't seem to figure out the translation to connect the array from the Image Orientation tag into a definitive statement regarding which plane the image is being taken from.

What I specifically am curious about is how I can take an array [0.5,0,-0.8660254,0,1,0] and be able to say it is a sagittal view like Roni says in this blog post: https://dicomiseasy.blogspot.com/2013/06/getting-oriented-using-image-plane.html especially when sagittal cuts are supposedly defined like this - ['0', '1', '0', '0', '0', '-1'].

I don't know if there is an existing plug in for Python users to be able to make this conversion, but if someone can help me understand how the process is completed, I'm happy to make one.

Mark Rotteveel
  • 100,966
  • 191
  • 140
  • 197
A.Code.1
  • 161
  • 1
  • 3
  • 10

2 Answers2

4

If you are only interested in the plane of the image, you basically need the direction of the image z axis. You get this by taking the cross product of the row and column image direction vectors (e.g. the ones that are given in ImageOrientationPatient). In the resulting vector you can define the main direction as the absolute maximum of the values. In your example you have:

>>> import numpy as np
>>> a = np.array([0.5,0,-0.8660254])
>>> b = np.array([0,1,0])
>>> np.cross(a, b)
array([ 0.8660254, -0.       ,  0.5      ])

The main direction in this case is the patient x axis (with a secondary direction being the patient z axis), meaning a sagittal view, with a secondary transverse direction (x -> sagittal, y -> coronal, z -> transverse/axial).

So to determine the main scan direction you could do something like this:

import numpy as np

ds = dcmread(filename)
image_ori = ds.ImageOrientationPatient
image_y = np.array([image_ori[0], image_ori[1], image_ori[2]])
image_x = np.array([image_ori[3], image_ori[4], image_ori[5]])
image_z = np.cross(image_x, image_y)
abs_image_z = abs(image_z)
main_index = list(abs_image_z).index(max(abs_image_z))
if main_index == 0:
    main_direction = "sagittal"
elif main_index == 1:
    main_direction = "coronal"
else:
    main_direction = "transverse"
print(main_direction)
MrBean Bremen
  • 14,916
  • 3
  • 26
  • 46
  • Thanks, this makes sense. I'm curious then about another example that Roni gave with an input that is [0.707, 0.707, 0, 0.707, 0.707, 0]. Taking the cross of [0.707, 0.707, 0] with [0.707, 0.707, 0] returns [0., 0., 0.] which the method will return as "sagittal" but Roni's post describes it as "transverse oblique" – A.Code.1 Jan 11 '22 at 21:13
  • I just had a look at that example, and I can't see how that could be right. I think this must be a mistake, the row and column vectors can't be the same. This should be something like [0.707, 0.707, 0, 0.707, - 0.707, 0] I think, e.g. he forgot a minus. Can have another look tomorrow when I get back to the PC. – MrBean Bremen Jan 11 '22 at 21:54
  • Actually check the comments in that post - they also point that out, and he agrees, just seems to have forgotten to adapt the post. – MrBean Bremen Jan 11 '22 at 21:58
  • Yep, the negative sign does resolve this, thank you! – A.Code.1 Jan 12 '22 at 00:16
1

MrBean has already answered this really well, but I thought I'd elaborate a bit more on your question.

You say: 'What I specifically am curious about is how I can take an array [0.5,0,-0.8660254,0,1,0] and be able to say it is a sagittal view'

The sagittal, axial, coronal views are, if you want to be strict about it, only those along the axes of the patient. So, as you point out, only those with all 1's or 0's in the IOP vectors.

So, for your example above, you *couldn't say [0.5,0,-0.8660254,0,1,0] is sagittal, as it's cross product isn't strictly in-line with any of the unit vectors.

But what everyone does, and what Roni in his great post somewhat says, is that you can say a particular orientation is more, or closest to axial, or sagittal or coronal, and call it that.

Or, you could put a tolerance around each of the unit views and only call it coronal or axial or sagittal if it's within that tolerance of the unit vectors - and otherwise say it's undefined. (you can think of this tolerance being like an angle - so, for example, if it's with 10 degrees of a unit vector).

So MrBean's code above is actually doing just that - for a given input IOP it's finding the closest unit vector to the normal (note the call to max):

main_index = list(abs_image_z).index(max(abs_image_z))

In my code, I use this definition regularly, so made my own method to calculate it. Note, like MrBean's code, it also will find the 'closest' orientation, iop_rounded = [round(x) for x in iop]. This will round all values to 0 or 1, and so at the end of the day, is equivalent to what MrBean used.

def get_orientation(dcm: DcmObj) -> str:
    """Get orientation of dicom object.

    Returns one of
        Transverse
        Coronal
        Sagittal
        NA  (if ImageOrientationPatient not available)
    """

    iop = getattr(dcm, 'ImageOrientationPatient', None)

    if iop is None:
        return 'NA'

    iop_rounded = [round(x) for x in iop]
    plane_cross = np.cross(iop_rounded[0:3], iop_rounded[3:6])
    plane = [abs(x) for x in plane_cross]

    if plane[0] == 1:
        return 'Sagittal'
    elif plane[1] == 1:
        return 'Coronal'
    elif plane[2] == 1:
        return 'Transverse'
    else:
        return 'NA'
Richard
  • 3,024
  • 2
  • 17
  • 40
  • Same question as I left on MrBean's response with regards to Roni giving an input that is [0.707, 0.707, 0, 0.707, 0.707, 0]. This method returns this as "NA" (cross product is all 0) but Roni's post describes it as "transverse oblique"? – A.Code.1 Jan 11 '22 at 21:16
  • OK - figured it out! I'd read his blog ages ago as well, and was thinking something was broken in my code!. The problem is he has his axes labelled incorrectly - they are missing a negative. In you above 2-vector, both vectors are the same: 0.7, 0.7, 0. That is how he has labelled them in his image. They are the same vector, so of course the cross product is zero. However this is incorrect (and if you look at the image you can see that) - one of the x or y numbers needs to be negative to make the vectors orthogonal. So try -0.7, 0.7, 0 and 0.7, 0.7, 0 - then it will be correct – Richard Jan 11 '22 at 23:26
  • Awesome, the negative sign resolves the error, thanks! As a small comment, I think by rounding before taking the cross product (iop_rounded before plane_cross), it sometimes throws off the calculation. For example with the -0.7, 0.7, 0 and 0.7, 0.7, 0, it returns NA as opposed to Transverse because the final array is [0, 0, 2]. Doing the rounding after the cross product seems to resolves this, but otherwise the solution is brilliant! – A.Code.1 Jan 12 '22 at 00:19
  • thanks for the catch! ...I only tested when the vectors are close to the unit vectors, so missed that. – Richard Jan 12 '22 at 00:33