# Get world coordinate of maximum intensity voxel

**URL:** <https://discourse.slicer.org/t/get-world-coordinate-of-maximum-intensity-voxel/4168>\
**Category:** Development\
**Created:** [September 22, 2018, 11:54am UTC](https://discourse.slicer.org/t/get-world-coordinate-of-maximum-intensity-voxel/4168 "2018-09-22T11:54:40Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![shahrokh](https://sea2.discourse-cdn.com/flex002/user_avatar/discourse.slicer.org/shahrokh/32/1499_2.png) [@shahrokh](https://discourse.slicer.org/u/shahrokh)\
**Post date:** [September 22, 2018, 11:54am UTC](https://discourse.slicer.org/t/get-world-coordinate-of-maximum-intensity-voxel/4168/1 "2018-09-22T11:54:41Z")

</div>

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](https://us1.discourse-cdn.com/flex002/uploads/slicer/original/3X/f/7/f75f780a6bb70811c419855bca601cd82c945ed2.png)

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

---

<div class="post-metadata">

**Author:** ![lassoan](https://sea2.discourse-cdn.com/flex002/user_avatar/discourse.slicer.org/lassoan/32/13_2.png) [@lassoan](https://discourse.slicer.org/u/lassoan)\
**Post date:** [September 22, 2018, 12:20pm UTC](https://discourse.slicer.org/t/get-world-coordinate-of-maximum-intensity-voxel/4168/2 "2018-09-22T12:20:05Z")

</div>

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](https://www.slicer.org/wiki/Documentation/Nightly/ScriptRepository#Get_volume_voxel_coordinates_from_markup_fiducial_RAS_coordinates), 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.

---

<div class="post-metadata">

**Author:** ![shahrokh](https://sea2.discourse-cdn.com/flex002/user_avatar/discourse.slicer.org/shahrokh/32/1499_2.png) [@shahrokh](https://discourse.slicer.org/u/shahrokh)\
**Post date:** [September 23, 2018, 2:27pm UTC](https://discourse.slicer.org/t/get-world-coordinate-of-maximum-intensity-voxel/4168/3 "2018-09-23T14:27:03Z")

</div>

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](https://us1.discourse-cdn.com/flex002/uploads/slicer/original/3X/5/7/5726ea77613f3d79e4b45521b93afeea38b5c5d9.png)

---

<div class="post-metadata">

**Author:** ![lassoan](https://sea2.discourse-cdn.com/flex002/user_avatar/discourse.slicer.org/lassoan/32/13_2.png) [@lassoan](https://discourse.slicer.org/u/lassoan)\
**Post date:** [September 23, 2018, 10:37pm UTC](https://discourse.slicer.org/t/get-world-coordinate-of-maximum-intensity-voxel/4168/4 "2018-09-23T22:37:57Z")

</div>

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.