Search code examples
pythondicompydicom

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.


Solution

  • 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)