Extrac Silhouette point clouds of each sagittal slice

Support
1

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
image940×607 98.7 KB

image
image1567×1062 67.4 KB

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.

2

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.

3

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?