Centroid determination

Support
1

I am trying to determine the centroid of several previously segmented materials. I can’t find this option in the Quantification module “Segment Statistics” of this software.
Is there any solution to this problem or the software 3D Slicer does’nt have any module that determines the centroid of an already segmented material?

I am currently using version 4.8.1 of the 3D Slicer software.

1 Like
2

I have the same issue as yours.

3

It is shown in this example how to get the centroid using numpy and then transform the raw array coordinates to scene RAS coordinates: https://www.slicer.org/wiki/Documentation/Nightly/ScriptRepository#Get_centroid_of_a_segment_in_world_.28RAS.29_coordinates

4

I have some follow-up questions:

  1. If there is only vtkMRMLModelNode, how can I find centroid of a VTK model?
  2. If 1) is not possible, how can I convert model to segmentation in a loadable script? I cannot find any associated function from slicer.modules.models.logic()

Thank you very much!

5

You can get center point very quickly as center of bounding box, with a little more computation from centroid of model points, or the exa t centroid (center of gravity assuming the model is filled with material of uniform density) using Segment statistics module.

6

I see, thank you.
I ended up implementing it without referring to bounding box, here’s my sample code as future reference for anyone who needs it:

md = slicer.util.getNode(‘ModelX’)
pd = md.GetPolyData().GetPoints().GetData()
from vtk.util.numpy_support import vtk_to_numpy
import numpy as np
print(np.average(vtk_to_numpy(pd), axis=0))

The output should print the average coordinate of all points within a vtkMRMLModelNode.

3 Likes
7

Thank you for sharing the solution! You can use slicer.util helper functions to simplify things a bit:

modelNode = slicer.util.getNode(‘ModelX’)
pointCoordinates = slicer.util.arrayFromModelPoints(modelNode)
import numpy as np
print(np.average(vtk_to_numpy(pd), axis=0))

or even shorter (a bit less explicit, so I would only recommend for testing and troubleshooting, not in modules):

import numpy as np
print(np.average(array('ModelX'), axis=0))
1 Like
8

Oh, great! I like the two line solution, and it reads quite explicit to me :slight_smile:

9

I have a question about getting a centroid of segment.
I ran this code below, but it couldn’t work.

segmentationNode = slicer.util.getNode("vtkMRMLSegmentationNode1")
segmentId = "Segment_1"

# Get array voxel coordinates
seg=slicer.util.arrayFromSegment(segmentationNode, segmentId)
# numpy array has voxel coordinates in reverse order (KJI instead of IJK)
# and the array is cropped to minimum size in the segmentation
mean_KjiCropped = [coords.mean() for coords in np.nonzero(seg)]

# Get segmentation voxel coordinates
segImage = segmentationNode.GetBinaryLabelmapRepresentation(segmentId)
segImageExtent = segImage.GetExtent()
# origin of the array in voxel coordinates is determined by the start extent
mean_Ijk = [mean_KjiCropped[2], mean_KjiCropped[1], mean_KjiCropped[0]] + np.array([segImageExtent[0], segImageExtent[2], segImageExtent[4]])

# Get segmentation physical coordinates
ijkToWorld = vtk.vtkMatrix4x4()
segImage.GetImageToWorldMatrix(ijkToWorld)
mean_World = [0, 0, 0, 1]
ijkToRas.MultiplyPoint(np.append(mean_Ijk,1.0), mean_World)
mean_World = mean_World[0:3]

# If segmentation node is transformed, apply that transform to get RAS coordinates
transformWorldToRas = vtk.vtkGeneralTransform()
slicer.vtkMRMLTransformNode.GetTransformBetweenNodes(segmentationNode.GetParentTransformNode(), None, transformWorldToRas)
mean_Ras = transformWorldToRas.TransformPoint(mean_World)

# Show mean position value and jump to it in all slice viewers
print(mean_Ras)
slicer.modules.markups.logic().JumpSlicesToLocation(mean_Ras[0], mean_Ras[1], mean_Ras[2], True)

from Documentation/Nightly/ScriptRepository - Slicer Wiki
I got error below

TypeError                                 Traceback (most recent call last)
Cell In[15], line 11
      8 mean_KjiCropped = [coords.mean() for coords in np.nonzero(seg)]
     10 # Get segmentation voxel coordinates
---> 11 segImage = segmentationNode.GetBinaryLabelmapRepresentation(segmentId)
     12 segImageExtent = segImage.GetExtent()
     13 # origin of the array in voxel coordinates is determined by the start extent

TypeError: GetBinaryLabelmapRepresentation() takes exactly 2 arguments (1 given)

Could you teach me the reason why this error happen?
Thank you

10

You should probably call the function as shown in the latest script repository.

1 Like