How to add Patient Position tag (Origin) to each slice in a SimpleITK image?

Beginner Questions
, , ,
1

I have a multi-slice volume [shape (256,256,14)] stored in a .vtk file. The header of the file, contains the following information:

# vtk DataFile Version 3.0    
vtk output
BINARY
DATASET STRUCTURED_POINTS
DIMENSIONS 256 256 14
SPACING 1.25 1.25 8
ORIGIN 111.81 -205.959 104.612
POINT_DATA 917504
SCALARS scalars float
LOOKUP_TABLE default

The ORIGIN refers to the Patient Position [DICOM Tag (0020,0032)] for slice 0.

Each slice [there are 14 slices], has a common Orientation Position [DICOM Tag (0020,0037)] but different Patient Position values.

I read this .vtk file as a SimpleITK image and use the TransformPhysicalPointToIndex() to convert world coordinates to image coordinates.

import SimpleITK as sitk

path = "multi_slice_volume.vtk"

reader = sitk.ImageFileReader()
reader.SetFileName(path)
multi_slice_vlm = reader.Execute()

world_coord = np.array([[72.424, -171.130, 125.821]])
image_coord = multi_slice_vlm.TransformPhysicalPointToIndex(world_coord)

I know in advance the correct value for image_coord - it should be (5,5,7) [this point belongs to slice 7]. However, the image_coord that I get is (13, 26, 4).

After some research, I realised this happens because, during the conversion world_coord->image_coord:

  • a) the Orientation Patient value used is the Simple ITK default value (1,0,0,0,1,0,0,0,1),
  • b) the Position Patient value used is the one from slice 0 (111.81 -205.959 104.612)

I was able to fix the a) issue using:

multi_slice_colume.SetDirection((0.43, -0.58, 0.0, 0.85, 0.014, 0.0, -0.29, -0.81, 1.0))

However, I’m struggling with issue b). For each slice of the Simple ITK image, I need to set the corresponding Position Patient tag. Here’s what I tried:

multi_slice_vlm[:, :, 0].SetOrigin([111.809579849243, -205.95926731824, 104.612316131591]) # slice 0
multi_slice_vlm[:, :, 1].SetOrigin( [106.316583633422, -201.75588673353, 108.632091522216]) # slice 1
    (...)
multi_slice_vlm[:, :, 13].SetOrigin( [100.823587417602, -197.55250614881, 112.651866912841]) # slice 13

But in the end, I continue to get the wrong value for image_coords.

2

You should use TransformIndexToPhysicalPoint() to get slice positions. So set proper direction to the image, then in a for loop construct indices [0,0,0] through [0,0,13]. For each index, get its physical position via multi_slice_colume.TransformIndexToPhysicalPoint(index).

3

Thank you @dzenanz. I’m quite new to ITK so I’m having some trouble understanding your answer.

Step 1 - “set proper direction to the image”:

path = "multi_slice_volume.vtk"

reader = sitk.ImageFileReader()
reader.SetFileName(path)
multi_slice_vlm = reader.Execute()

multi_slice_vlm.SetDirection((0.43, -0.58, 0.0, 0.85, 0.014, 0.0, -0.29, -0.81, 1.0))

Step 2 - " in a for loop construct indices [0,0,0] through [0,0,13]"
I checked the link you provided (namely the Python example) but I was not able to understand how to use the itk.Index().

Also, it seems that I need to first convert the SimpleITK image in a ITK image, right?