Load DICOM Surface Segmentation Objects

Support Feature requests
1

Slicer version: 4.6.2

Hi,
I just wondered, if there is any possibility to load DICOM Surface Segmentation Objects, as described in the DICOM Supplement 132 with the 3D Slicer.

Thank you in advance! :slight_smile:
Best regards
Julia

2

Hi Julia -

I wrote a first version that worked well on a sample dataset:

it should work to add this as loadable module.

I didn’t have any sharable sample data at the time so I didn’t add this to the core. But I’d like to finish this sometime so if you have any example data to share testing/debugging that would be great.

3

I did not yet get the plugin working in slicer but managed to produce working code with the massive help of the author of the aliza dicom viewer.

It imports a dicom surface segmentation and writes it to a vtk and stl file. The key part was to subtract 1 for the indices

polys.InsertCellPoint(index-1)

The dicom file to test this was created on a Siemens syngo workstation. I hope this helps. Would love to see this working in Slicer.

import glob, os, json
from datetime import datetime
import string
import vtk
import logging


def load(loadable):
    import dicom
    dataset = dicom.read_file(loadable)
    loadSurfaceSegmentationDataset(dataset)

def loadSurfaceSegmentationDataset(dataset):
    for surface in dataset.SurfaceSequence:
      coordinates = surface.SurfacePointsSequence[0].PointCoordinatesData

      # generate points
      points = vtk.vtkPoints()
      polyData = vtk.vtkPolyData()
      polyData.SetPoints(points)

      for pointIndex in xrange(len(coordinates)/3):
        pointStart = pointIndex*3
        point = coordinates[pointStart:pointStart+3]
        point[0] *= -1
        point[1] *= -1
        points.InsertNextPoint(*point)

      polys = vtk.vtkCellArray()
      polyData.SetPolys(polys)

      # generate triangles
      extremes = [100,1]
      triIndices = surface.SurfaceMeshPrimitivesSequence[0].TrianglePointIndexList
      triangles = len(triIndices)/6
      for triangle in range(triangles):
        polys.InsertNextCell(3)
        for vertex in range(3):
          vertexIndex = 6*triangle + 2*vertex;
          low = ord(triIndices[vertexIndex])
          high = ord(triIndices[vertexIndex+1])
          index = (high << 8) + low
          polys.InsertCellPoint(index-1)

          if index < extremes[0]:
            extremes[0] = index
          if index > extremes[1]:
            extremes[1] = index

    # visualization
    mapper = vtk.vtkPolyDataMapper()
    mapper.SetInputData(polyData)

    print polyData

    actor = vtk.vtkActor()
    actor.SetMapper(mapper)
    actor.GetProperty().SetColor(255,255,255)

    ren=vtk.vtkRenderer()
    renWin=vtk.vtkRenderWindow()
    renWin.AddRenderer(ren)
    iren=vtk.vtkRenderWindowInteractor()
    iren.SetRenderWindow(renWin)

    ren.AddActor(actor)
    ren.SetBackground(0.1, 0.2, 0.4)
    renWin.SetSize(200, 200)

    iren.Initialize()

    ren.ResetCamera()
    ren.ResetCameraClippingRange()

    renWin.SetSize(800,600)
    renWin.Render()

    iren.Start()

    # build vtk file
    writer = vtk.vtkPolyDataWriter()
    writer.SetInputData(polyData)
    output_path = "test"
    writer.SetFileName(output_path + '.vtk')
    writer.Write()

    # build stl file
    surface_filter = vtk.vtkDataSetSurfaceFilter()
    surface_filter.SetInputData(polyData)

    triangle_filter = vtk.vtkTriangleFilter()
    triangle_filter.SetInputConnection(surface_filter.GetOutputPort())

    writer = vtk.vtkSTLWriter()
    writer.SetFileName(output_path + '.stl')
    writer.SetInputConnection(triangle_filter.GetOutputPort())
    writer.SetFileTypeToBinary()
    writer.Write()

###########################################################################################
if __name__ == "__main__":

    load('LASegSample.dcm')
4

Just one more thing. There is a surfaces.dcm out in the wild that serves as as test file. It is referenced in this thread.

https://groups.google.com/forum/#!topic/comp.protocols.dicom/3mhBWlLZSuM

Be careful with this file. Indices seem to be off. To avoid confusion always work with files directly exported from software such as Syngo (Siemens).

5

Hi.
Confirm. Indices must start with 1, not with 0 like in ‘surfaces.dcm’ i found somewhere and first used as reference too, Problem is too few example files. Many thanks to Sebastian Hilbert for sharing a real sample from Siemens Workstation.

DICOM PS3.3 C.27.2.1.1 Point Coordinates Data
“When referencing individual points the index of the first point shall be 1.”
http://dicom.nema.org/medical/dicom/current/output/chtml/part03/sect_C.27.2.html#sect_C.27.2.1.1

And ‘surfaces.dcm’ hat 32 bit int in Triangle Point Index List OW VR too. Both.

BTW. Attribute Triangle Point Index List is retired.
PS3.3 C.27.4.1 Surface Mesh Primitives Macro Attribute Descriptions
“In a previous edition, other Attributes were used that had an OW VR and a limitation to no more than 65535 points per surface. These have been retired and replaced with new Attributes. See PS 3.3 2014a.”
New attribute for triangle index list is “Long Triangle Point Index List” (0066,0041). But is not so critical.

Best Regards,
M

1 Like
6

Siemens SEG LASegSample.dcm

1 Like