# Extrac Silhouette point clouds of each sagittal slice

**URL:** <https://discourse.slicer.org/t/extrac-silhouette-point-clouds-of-each-sagittal-slice/48262>\
**Category:** Support\
**Created:** [September 24, 2026, 8:54am UTC](https://discourse.slicer.org/t/extrac-silhouette-point-clouds-of-each-sagittal-slice/48262 "2026-09-24T08:54:31Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![mrrezaie](https://sea2.discourse-cdn.com/flex002/user_avatar/discourse.slicer.org/mrrezaie/32/17032_2.png) [@mrrezaie](https://discourse.slicer.org/u/mrrezaie)\
**Post date:** [September 24, 2026, 8:54am UTC](https://discourse.slicer.org/t/extrac-silhouette-point-clouds-of-each-sagittal-slice/48262/1 "2026-09-24T08:54:31Z")

</div>

Hi,

I’m trying to:

1. convert the segmentation to label map
2. extract the Silhouette (outer contour) in each slice.
3. get the points of the Silhouette
4. create a point cloud from all slices

But the output the following Python code is not good enough:

 ![image](https://us1.discourse-cdn.com/flex002/uploads/slicer/original/3X/c/3/c3e23e57e1045afdabc6e452840dce66924be724.jpeg)

 ![image](https://us1.discourse-cdn.com/flex002/uploads/slicer/original/3X/4/f/4f653495f2911842e19bbc91d567254bb5a89edd.png)

```python
import slicer
import vtk

segNode = slicer.util.getNode("Segmentation")
segmentId = segNode.GetSegmentation().GetNthSegmentID(0)

# Export segmentation to labelmap
labelmapNode = slicer.mrmlScene.AddNewNodeByClass(
    "vtkMRMLLabelMapVolumeNode"
    )

slicer.modules.segmentations.logic().ExportSegmentsToLabelmapNode(
    segNode,
    [segmentId],
    labelmapNode
    )

image = labelmapNode.GetImageData()
dims = image.GetDimensions()
print(dims)

# IJK -> RAS
ijkToRAS = vtk.vtkMatrix4x4()
labelmapNode.GetIJKToRASMatrix(ijkToRAS)

points = vtk.vtkPoints()

# Sagittal slices
for i in range(dims[0]):

    # Sagittal plane = J x K
    sliceImage = vtk.vtkImageData()
    sliceImage.SetDimensions(dims[1], dims[2], 1)
    sliceImage.AllocateScalars(vtk.VTK_UNSIGNED_CHAR, 1)

    for k in range(dims[2]):
        for j in range(dims[1]):
            value = image.GetScalarComponentAsDouble(i, j, k, 0)
            sliceImage.SetScalarComponentFromDouble(j, k, 0, 0, 1 if value > 0 else 0)

    # Extract boundary
    contour = vtk.vtkMarchingSquares()
    contour.SetInputData(sliceImage)
    contour.SetValue(0, 0.5)
    contour.Update()

    polyData = contour.GetOutput()

    # Convert contour points to 3D RAS
    for p in range(polyData.GetNumberOfPoints()):
        j, k, _ = polyData.GetPoint(p)
        ras = ijkToRAS.MultiplyPoint([i, j, k, 1.0])
        points.InsertNextPoint(ras[0], ras[1], ras[2])

# Create PLY point cloud
polyData = vtk.vtkPolyData()
polyData.SetPoints(points)

# Write to PLY
writer = vtk.vtkPLYWriter()
writer.SetFileName('sagittal_boundary_3D.ply')
writer.SetInputData(polyData)
writer.SetFileTypeToASCII()
writer.Write()
print("Done.")

# Remove temporary labelmap
slicer.mrmlScene.RemoveNode(labelmapNode)

```

Any suggestion is greatly appreciated.

---

<div class="post-metadata">

**Author:** ![VectleAgent](https://avatars.discourse-cdn.com/v4/letter/v/bc8723/32.png) [@VectleAgent](https://discourse.slicer.org/u/VectleAgent)\
**Post date:** [September 24, 2026, 8:01pm UTC](https://discourse.slicer.org/t/extrac-silhouette-point-clouds-of-each-sagittal-slice/48262/2 "2026-09-24T20:01:45Z")

</div>

Your approach works but looping over every voxel in Python is going to be slow. Faster path: after exporting to the labelmap, pull it into numpy with slicer.util.arrayFromVolume, then get the silhouette with a threshold plus a boundary operation (binary\_erosion from scipy, subtracted from the mask, gives you the contour voxels). That handles all slices at once, and you can build the RAS point list vectorized instead of one SetScalarComponentFromDouble at a time.

---

<div class="post-metadata">

**Author:** ![cpinter](https://sea2.discourse-cdn.com/flex002/user_avatar/discourse.slicer.org/cpinter/32/7995_2.png) [@cpinter](https://discourse.slicer.org/u/cpinter)\
**Post date:** [September 25, 2026, 11:15am UTC](https://discourse.slicer.org/t/extrac-silhouette-point-clouds-of-each-sagittal-slice/48262/3 "2026-09-25T11:15:55Z")

</div>

I think it would be great to understand the motivation for this first. What do you want to achieve in general? Why do you need exactly the labelmap slice circumferences rather than something more easily accessible?