OrientImageFilter works differently before and after v4.13.1.post1?

[sbonaretti]

Hi!

I have a function to orient images to RAI using OrientImageFilter (Please, find images and whole function below). The same function works if I use ITK 4.13.1.post1, but it does not work if I use ITK versions after that. Not working means that the image does not change orientation and I do not get any error message, like nothing happens. Do you have any suggestions on what could cause the difference?

This is actually an issue because I cannot update to new ITK versions for my tool (pyKNEEr).

FYI: I am also part of an issue about this in SimpleITK because my whole package is written in SimpleITK, except this function. It seems like the reason is incompatibility issues between dicom and ITK.

Thanks a lot!
Serena

def orientation_to_rai(img):

    # this function is implemented using ITK since SimpleITK has not implemented the filter "OrientImageFilter" yet
    # the ITK orientation system is from this blog post: https://itk.org/pipermail/insight-users/2017-May/054606.html
    # comparison ITK - simpleITK filters: https://itk.org/SimpleITKDoxygen/html/Filter_Coverage.html
    # see also: https://github.com/fedorov/lidc-idri-conversion/blob/master/seg/seg_converter.py

    # change image name to specify it is the simpleITK one
    img_sitk = img

    # get characteristics of simpleITK image
    size_in_sitk      = img_sitk.GetSize()
    spacing_in_sitk   = img_sitk.GetSpacing()
    origin_in_sitk    = img_sitk.GetOrigin()
    direction_in_sitk = img_sitk.GetDirection()

    # allocate ITK image (type and size)
    Dimension   = 3
    PixelType   = itk.F
    ImageTypeIn = itk.Image[PixelType, Dimension]
    img_itk = ImageTypeIn.New()
    sizeIn_itk = itk.Size[Dimension]()
    for i in range (0,Dimension):
        sizeIn_itk[i] = size_in_sitk[i]
    region = itk.ImageRegion[Dimension]()
    region.SetSize(sizeIn_itk)
    img_itk.SetRegions(region)
    img_itk.Allocate()

    # pass image from simpleITK to numpy
    img_py  = sitk.GetArrayFromImage(img_sitk)

    # pass image from numpy to ITK
    img_itk = itk.GetImageViewFromArray(img_py)

    # pass characteristics from simpleITK image to ITK image (except size, assigned in allocation)
    spacing_in_itk = itk.Vector[itk.F, Dimension]()
    for i in range (0,Dimension):
        spacing_in_itk[i] = spacing_in_sitk[i]
    img_itk.SetSpacing(spacing_in_itk)

    origin_in_itk  = itk.Point[itk.F, Dimension]()
    for i in range (0,Dimension):
        origin_in_itk[i]  = origin_in_sitk[i]
    img_itk.SetOrigin(origin_in_itk)

    direction_in_itk = itk.Matrix[itk.F,3,3]()
    direction_in_itk = img_itk.GetDirection().GetVnlMatrix().set(0,0,direction_in_sitk[0]) # r,c,value
    direction_in_itk = img_itk.GetDirection().GetVnlMatrix().set(0,1,direction_in_sitk[1])
    direction_in_itk = img_itk.GetDirection().GetVnlMatrix().set(0,2,direction_in_sitk[2])
    direction_in_itk = img_itk.GetDirection().GetVnlMatrix().set(1,0,direction_in_sitk[3]) # r,c,value
    direction_in_itk = img_itk.GetDirection().GetVnlMatrix().set(1,1,direction_in_sitk[4])
    direction_in_itk = img_itk.GetDirection().GetVnlMatrix().set(1,2,direction_in_sitk[5])
    direction_in_itk = img_itk.GetDirection().GetVnlMatrix().set(2,0,direction_in_sitk[6]) # r,c,value
    direction_in_itk = img_itk.GetDirection().GetVnlMatrix().set(2,1,direction_in_sitk[7])
    direction_in_itk = img_itk.GetDirection().GetVnlMatrix().set(2,2,direction_in_sitk[8])

    # make sure image is float for the orientation filter (GetImageViewFromArray sets it to unsigned char)
    ImageTypeIn_afterPy = type(img_itk)
    ImageTypeOut        = itk.Image[itk.F, 3]
    CastFilterType      = itk.CastImageFilter[ImageTypeIn_afterPy, ImageTypeOut]
    castFilter          = CastFilterType.New()
    castFilter.SetInput(img_itk)
    castFilter.Update()
    img_itk             = castFilter.GetOutput()

    # define ITK orientation system  (from the blog post: https://itk.org/pipermail/insight-users/2017-May/054606.html)
    ITK_COORDINATE_UNKNOWN   = 0
    ITK_COORDINATE_Right     = 2
    ITK_COORDINATE_Left      = 3
    ITK_COORDINATE_Posterior = 4
    ITK_COORDINATE_Anterior  = 5
    ITK_COORDINATE_Inferior  = 8
    ITK_COORDINATE_Superior  = 9
    ITK_COORDINATE_PrimaryMinor   = 0
    ITK_COORDINATE_SecondaryMinor = 8
    ITK_COORDINATE_TertiaryMinor  = 16
    ITK_COORDINATE_ORIENTATION_RAI = ( ITK_COORDINATE_Right    << ITK_COORDINATE_PrimaryMinor ) \
                                   + ( ITK_COORDINATE_Anterior << ITK_COORDINATE_SecondaryMinor ) \
                                   + ( ITK_COORDINATE_Inferior << ITK_COORDINATE_TertiaryMinor )

    # change orientation to RAI
    OrientType = itk.OrientImageFilter[ImageTypeOut,ImageTypeOut]
    filter     = OrientType.New()
    filter.UseImageDirectionOn()
    filter.SetDesiredCoordinateOrientation(ITK_COORDINATE_ORIENTATION_RAI)
    filter.SetInput(img_itk)
    filter.Update()
    img_itk    = filter.GetOutput()


    # get characteristics of ITK image
    spacing_out_itk   = img_itk.GetSpacing()
    origin_out_itk    = img_itk.GetOrigin()

    # pass image from itk to numpy
    img_py = itk.GetArrayViewFromImage(img_itk)

    # pass image from numpy to simpleitk
    img = sitk.GetImageFromArray(img_py)

    # pass characteristics from ITK image to simpleITK image (except size, implicitely passed)
    spacing = []
    for i in range (0, Dimension):
        spacing.append(spacing_out_itk[i])
    img.SetSpacing(spacing)

    origin = []
    for i in range (0, Dimension):
        origin.append(origin_out_itk[i])
    img.SetOrigin(origin)

    direction = []
    direction.append(1.0)
    direction.append(0.0)
    direction.append(0.0)
    direction.append(0.0)
    direction.append(1.0)
    direction.append(0.0)
    direction.append(0.0)
    direction.append(0.0)
    direction.append(1.0)
    img.SetDirection(direction)

    return img

[dzenanz]

No significant changes were done to this class since be3c980 on 2016-06-19, which was before 4.13.1 tag. Change was caused by something else, perhaps DICOM reading? In your debugging Jupyter notebook, can you also print image's direction?

[blowekamp]

We have been struggling the best way to support this filter in SimpleITK. I also reported the discrepancy of the orientation here: #1163 Do to compatibility, I can't for see the current itk::OrientImageFilter changing it's definitions to match DICOM.

When the OrientImageFilter is used, frequently the ResampleImageFilter may be a better choice, because the OrientImageFilter does not handle matrices with rotation. So you code could be made more robust by changing to use the ResampleImageFilter and correctly specifying the desired meta data and the change in the output geometry of the image.

The solution to the OrientImageFilter is to create a new DICOMOrientImageFilter, with the correct
orientations which match the DICOM standard. I'll am going to work on creating this new filter in the itkSimpleITKFilters remote module today.

[sbonaretti]

I have dug deeper into my code and I have found that the issue is in the allocation of the ITK image, and specifically in setting the direction!
Because of the way my code is structured, I have to read the image in SimpleITK. So, to use the ITK filter, I move it to Numpy, then to ITK. In these chain, I manually assign the direction of the SimpleITK image to the ITK image. I have just found here that there is a new way to set direction in ITK:

img.SetDirection(itk.matrix_from_array(3x3_numpy_array))

That fixed the problem! My ITK image was keeping the default direction instead of having the one that (I thought) I was assigning.

Thanks a lot for your kind replies! And thanks a lot for your effort to make the filter available in SimpleITK!!!
Serena