Get center of gravity of a segment

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.

Here you go:

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 ?

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.

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)

See this post: How to display and move the origin coordinates of the surface model (.STL)

1 Like

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

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()

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

Thanks for explaining!

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

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.

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

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

Through the UI, you can import the labelmap as segmentation, turn them into individual segments using island tool, and then use the segment statistics module to get the centroids for each of those segments.

If you will do it more than once, scripting is probably the better approach.

1 Like

Hi muratmaga,
I got the results like this

How can I get the position coordinates for each of the centroide of a segment and then replace it with a fiducial

thank you

regards,
Saima safdar

Hi,
Thanks I found the solution

You can find a complete example of placing markups at segment centroids: https://www.slicer.org/wiki/Documentation/Nightly/ScriptRepository#Get_centroid_of_each_segment

Here is another example, which draws oriented bounding boxes around each segment: https://www.slicer.org/wiki/Documentation/Nightly/ScriptRepository#Get_size.2C_position.2C_and_orientation_of_each_segment

1 Like

Hi Andras,
The electrodes which I extracted now I can represent them with fiducials.

is there a way to project these fiducials onto the brain cortical surface on the T1 MRI? I need to calculate the distance between the projected and actual electrode positions in CT scan.

Any idea how can I get? Can I use the nearest neighbour to find cells on brain and then get cell centroid and convert to fiducials or is there a better way of projecting onto brain surface.

Thank you

Regards,
Saima Safdar

Projection to a curved surface is not a well-defined operation: it can be done many different ways, with slightly different results.

Can you post a few screenshots that show where the electrode positions are, how the cortical surface looks like, and where do you think they should be projected to?