CODEWITHSUNDEEP
pythondicompydicom
A.Code.1
A.Code.1

Reputation: 161

Translate Image Orientation into axial, sagittal, or coronal plane

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.

Upvotes: 0

Views: 4027

Answers (2)

Richard
Richard

Reputation: 3394

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'

Upvotes: 1

MrBean Bremen
MrBean Bremen

Reputation: 16815

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)

Upvotes: 4

Related Questions