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.

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?