Get world coordinate of maximum intensity voxel

Development
1

Hi
Dear Developers and Users
I want to create sphere model centerd in the voxel with maximum value.

Steps:
1- Load PET images using “DICOM Browser”.
2- Click on “Center Volume” in the module of “Volumes”.
3- Enter the following commands:
PETnode = slicer.util.getNode(‘PETnode’)
PETnumpy = slicer.util.arrayFromVolume(PETnode)
import numpy as np
indices = np.where(PETnumpy == PETnumpy.max())
sphere = vtk.vtkSphereSource()
sphere.SetRadius(10.)
sphere.SetCenter(indices)
sphere.SetPhiResolution(30)
sphere.SetThetaResolution(30)
modelsLogic = slicer.modules.models.logic()
model = modelsLogic.AddModel(sphere.GetOutput())
model.GetDisplayNode().SetSliceIntersectionVisibility(True)
model.GetDisplayNode().SetSliceIntersectionThickness(1)
model.GetDisplayNode().SetColor(1,0,0)
sphere.Update()

As mentioned in the following screen shot, the position of the sphere is not correct (out of body).

Screenshot%20from%202018-09-22%2016-19-20

Please guide me. Where do I wrong it?
Please guide me.
Shahrokh

2

Numpy returns (K, J, I) coordinates (voxel I, J, K coordinates in reverse order), while fiducials are defined in anatomical (R, A, S) coordinates.

(R, A, S) to (I, J, K) conversion is described here, you need to perform the inverse of this operation. If you implemented the inverse script then please copy it here and we’ll post it to the script repository for future reference.

3

Dear Andras

Thanks a lot for your guidance. At now, I can center my sphere model in the voxel with maximum uptake. As mentioned you, firstly I changed indices array (indicesKJI, coordinates position of maximum uptake) from (K J I) to (I J K). After doing it, I calculate inverse matrix volumeRasToIjk with the name of volumeIjkToRas. Finally I applied volumeIjkToRas to indicesIJK.

Code written for doing it:

#Load PET images with "DICOM" module
#Rename PET node in "Data" to PETnode
#Click on "Center Volume" in "Volumes" node
PETnode = slicer.util.getNode('PETnode')
PETnumpy = slicer.util.arrayFromVolume(PETnode)
#PETnumpy.max()
import numpy as np
indicesKJI = np.where(PETnumpy == PETnumpy.max())  
indicesIJK = [indicesKJI[2], indicesKJI[1], indicesKJI[0]]
indicesIJK = np.append(indicesIJK,1)
volumeNode = getNode('PETnode')
volumeRasToIjk = vtk.vtkMatrix4x4()
volumeNode.GetRASToIJKMatrix(volumeRasToIjk)
import numpy as np
from numpy.linalg import inv
volumeRasToIjk = np.array([[volumeRasToIjk.GetElement(0,0), volumeRasToIjk.GetElement(0,1), volumeRasToIjk.GetElement(0,2), volumeRasToIjk.GetElement(0,3)],[volumeRasToIjk.GetElement(1,0), volumeRasToIjk.GetElement(1,1), volumeRasToIjk.GetElement(1,2), volumeRasToIjk.GetElement(1,3)],[volumeRasToIjk.GetElement(2,0), volumeRasToIjk.GetElement(2,1), volumeRasToIjk.GetElement(2,2), volumeRasToIjk.GetElement(2,3)],[volumeRasToIjk.GetElement(3,0), volumeRasToIjk.GetElement(3,1), volumeRasToIjk.GetElement(3,2), volumeRasToIjk.GetElement(3,3)]])
volumeIjkToRas = inv(volumeRasToIjk)
indicesRAS = np.dot(volumeIjkToRas, indicesIJK)
sphere = vtk.vtkSphereSource()
sphere.SetRadius(10.)
sphere.SetCenter(indicesRAS[0], indicesRAS[1], indicesRAS[2])
sphere.SetPhiResolution(30)
sphere.SetThetaResolution(30)
modelsLogic = slicer.modules.models.logic()
model = modelsLogic.AddModel(sphere.GetOutput())
model.GetDisplayNode().SetSliceIntersectionVisibility(True)
model.GetDisplayNode().SetSliceIntersectionThickness(3)
model.GetDisplayNode().SetColor(1,0,0)
sphere.Update()

Screenshot result:
Screenshot%20from%202018-09-23%2018-53-35

2 Likes
4

Thanks for sharing your results.

You may consider adding markup fiducial points at these positions. You can set the markup shape to sphere and make it as large as you prefer, and you can easily jump to the point positions in Markups module by enabling “Click to jump slices” and clicking on a markup in the list.