Dear All,
I have a vtkPolyData file (extension.vtk) with 12 points:
These points are the annotation of a 3D image data (more exactly, a multi-slice DICOM image - one folder with 14 slices/frames .dcm).
I need to convert each of the previous points to image coordinates. I.e. given a point (e.g. -177.747, -136.373, 21.174) I need to find the correspondent (i, j, k) image coordinates (where k is the slice index).
This is what I tried:
# load 3D volume
volume_reader = sitk.ImageSeriesReader()
dicom_names = volume_reader.GetGDCMSeriesFileNames(multi_slice_volume_path)
volume_reader.SetFileNames(dicom_names)
image = volume_reader.Execute()
# load vtk points
reader = vtk.vtkPolyDataReader()
reader.SetFileName(lmks_vtk_file_path)
reader.Update()
vtkdata = reader.GetOutput()
vtk_points = numpy_support.vtk_to_numpy(vtkdata.GetPoints().GetData())
# convert vtk points to image coordinates
image_points = np.zeros((vtk_points.shape))
for i, coord in enumerate(vtk_points):
coord = (float(coord[0]), float(coord[1]), float(coord[2]))
image_points[i] = np.array(image.TransformPhysicalPointToIndex(coord)).astype(int)
But the image coordinates that I get are not correct.
[[-33. 190. 24.]
[-52. 194. 23.]
[-45. 178. 21.]
[-35. 143. 26.]
[-42. 166. 29.]
[-72. 168. 26.]
[-60. 147. 24.]
[-25. 161. 24.]
[-35. 181. 27.]
[-66. 180. 25.]
[-53. 161. 23.]
[-29. 179. 23.]]
I get negative i coordinates (points outside the image boundaries) and the k coordinates are not within the expected range 0-13 (since I have a total of 14 slices).
I did further research but I couln’t find the solution. The closest clue was this post https://discourse.vtk.org/t/mapping-between-coordinate-position-x-y-z-and-voxel-index/5357/4 by @dgobbi
Any help will be most appreciated. Thanks.