Get center of gravity of a segment

Development
1

Hi,

I’m new on 3Dslicer with a new project, I try to search on topics but I’m still in troubles.

I’m using segmentation from CT-scan to isolate 2 masses.
I need you to know how I can have the gravity mass center. Can I have coordinates x/y/z ?
Does it run with models ? volumes ?

Thanks for all

Sorry for the mistakes I don’t speak english very well

2

use this code or something similar:

   n = getNode('your node') 
   s = n.GetSegmentation()
  ss = s.GetSegment('Segment_7')
  pd = ss.GetRepresentation('Closed surface')
  com = vtk.vtkCenterOfMass()
  com.SetInputData(pd)
  com.Update()
  com.GetCenter()
3

Thanks !

The problem I have is that :

With this line pd = ss.GetRepresentation(‘Closed surface’) it requires a Closed surface representation but when I segment my 2 masses with threshold (and split island) it works with binary labelmap.

How can I do ?

(So it’s good with painted segments, but not tresholded)

Thanks

4

You can also use segmentation node’s GetSegmentCenter or GetSegmentCenterRAS method. It works for both labelmap and closed surface representations. It is used for giving a position inside the segment (used when you right-click the segment and select “Jump slices”) and so it does not use center of gravity. If you need strictly center of gravity then you can have a look at the source code of this method and create a modified version.

5

it works with binary labelmap.

You can convert from labelmap to segmentation node and vice versa using segmentations module.

6

If you have trouble modifying the GetSegmentCenter function to get the center of mass, this might help. In this Slicer programming tutorial the center of mass of a labelmap is calculated in python:
https://github.com/PerkLab/PerkLabBootcamp/blob/master/Doc/day3_2_SlicerProgramming.pptx

1 Like
7

It just occurred to me that we changed the bootcamp task to center of mass of fiducials not labelmap. Here’s an older one with the labelmap code:
https://www.dropbox.com/s/ee58jil89rw0qxd/20160504_SlicerProgramming_Pinter.pptx?dl=0 (slide 28). I need to reimplement just this for a completely unrelated reason, so if you want a few hours I’ll link a C++ code that calculates center of mass for a segment’s labelmap.

8

Here you go:

9

Thanks to everybody !
It’s now ok for me for center of gravity with your help !

Another question :

–> Can I analyse surface points coordinates to compute the most distants points for visualizing some 2 or 3 largest diameters of the mass ?

10

It would be more robust to compute principal axes from labelmap representation. You can then get diameters along those axes.

There may be ITK or VTK filters that can help. Let us know about your progress, we may give further help and we could integrate the final solution to Segment Statistics module.

11

Related question:

how can I change the origin of a model to the center of mass?

I found two useful functions GetBounds and GetMesh().GetCenter() . I think the origin is stored in the minimum values Xmin,Ymin and Zmin. The problem I can not find something like model.SetBounds. Here is the idea:

                 bounds =[0]*6
                 model.GetBounds(bounds)
                 center = model.GetMesh().GetCenter() 
                 bounds [0] = center[0] 
                 bounds [2 ]= center[1]
                 bounds [4] = center[2]
                #model.SetBounds(bounds)
12

See this post: How to display and move the origin coordinates of the surface model (.STL) - #2 by pieper

1 Like
13

I just found today that values are different when computing the center from a model e.g.:

 center = model.GetMesh().GetCenter()
 output: (-4.173086166381836, -155.35969924926758, -259.25328826904297)

and from computing the center using segment

 center = vtk.vtkCenterOfMass()
 center.SetInputData(pd)
 center.Update()
 center.GetCenter()
 output: (-4.462736238549897, -146.88848545865304, -257.3132177804882)

Apparently, the second one is more accurate as I get the same value in Blender. Is there a way to get the same values using model?

1 Like
14

Solved!

It seems I can get polydata from the model and use the second method e.g.

pd= model.GetPolyData()
center = vtk.vtkCenterOfMass()
center.SetInputData(pd)
center.Update()
center.GetCenter()
15

You have many options to compute “center” of a segment.

  1. Model node’s GetMesh().GetCenter() method returns center of the bounding box. Depending on the model’s shape it may be close to the center of mass, but they are computed very differently.
See method help 
>>> help(modelNode.GetMesh().GetCenter)
Help on built-in function GetCenter:

GetCenter(...)
    V.GetCenter() -> (float, float, float)
    C++: double *GetCenter()
    V.GetCenter([float, float, float])
    C++: void GetCenter(double center[3])
    
    Get the center of the bounding box. THIS METHOD IS NOT THREAD
    SAFE.

  1. vtk.vtkCenterOfMass() applied to a segment’s closed surface representation returns center of mass of the surface points. Depending on the shape and density of points on the surface, it may be similar to center of mass of the of the solid object enclosed by the surface, but they are not the same.

  2. You can compute true center of mass of the enclosed volume from binary labelmap representation (as @cpinter described above).

3 Likes
16

Thanks for explaining!

17

Thanks for asking this questions and for the answers, it helped me a lot. However, here’s a little note, in case someone else gets stuck like I did just now. I tried this code segment and this line:

pd = ss.GetRepresentation(‘Closed surface’)

returned a “None” for me. In this case, press “Show 3D” in the Segment Editor at least once, then the “Closed surface” representation gets calculated internally and this line returns a vtkPolyData object (a surface model) as expected. Alternatively, the calculation of this representation can also be triggered in code, via the following line (beware, I did not fully test whether the result makes sense and can be worked with further, I just checked whether pd is not None in the next line):

ss.AddRepresentation(‘Closed surface’,vtk.vtkPolyData())

Another tip: the following line assumes that the segment ID is known:

ss = s.GetSegment(‘Segment_7’)

If this ID is not known, the following line gets the segment by its Name, rather than by its ID:

ss = s.GetSegment(s.GetSegmentIdBySegmentName(‘Calcified mass’))

1 Like
18

Please do not do this, it will result in inconsistent content in the segmentation object, and subsequent operations may yield unwanted results. This is how you can trigger conversion internally:

s.CreateRepresentation(‘Closed surface’)

If you want to learn more about the infrastructure, then you can refer to the documentation:
https://www.slicer.org/wiki/Documentation/Nightly/Modules/Segmentations
or the source code, which is well commented.

19

Or you can use segmentationNode.CreateClosedSurfaceRepresentation() convenience function.

20

Hi Andras,
How can I get center of gravity or centriod position for each of these electrode positions. The dots are in one labelmap segmentation.

image

Thank you

Regards,
Saima Safdar