View DICOM SEG mesh with VTK and pydicom

loadSEGmesh.py

# View DICOM SEG Surface Segmentation.
#
# Run standalone (VTK Python 3 and pydicom >= 2 are required)
# python3 loadSEGmesh.py <file>
#
# or from Slicer >= 4.11
# Slicer --python-script loadSEGmesh.py <file>
# (caution with Slicer RAS space, s. options below)
#
#
# Example file (original Siemens):
# https://drive.google.com/file/d/1SmkP93YHf0Sq1pYBNjFJyewmAhdP1Omk/view?usp=sharing
#
#
#
# Options:
#
#
# Flip X, Y axis (Slicer specific)
dicom_lps_to_slicer_ras = False
#
# Write
write_vtk = False
write_stl = False
output_prefix = "surface-"
#
# Indices must start with 1,
# but bad datasets exist with first index 0.
index_start_with_0 = False
#
# Indices in TrianglePointIndexList (VR OW) must be 16 bits words (uint16),
# but bad datasets exist with 32 bit indices.
# Note: indices in LongTrianglePointIndexList (VR OL) are 32 bits dwords (uint32),
# the option has no effect at LongTrianglePointIndexList.
force_index_type_int = False
#
#
#
#
#
#
#
#
import glob
import os
import sys
import string
import random
import vtk
from struct import unpack

# CIELab to RGB
# https://getreuer.info/posts/colorspace/
def lab_inverse(x):
    if x >= 0.206896551724137931:
        y = x * x * x
    else:
        y = (108.0 / 841.0) * (x - (4.0 / 29.0))
    return y

def gamma_correction(x):
    if x <= 0.0031306684425005883:
        y = 12.92 * x
    else:
        y = 1.055 * pow(x, 0.416666666666666667) - 0.055
    return y

# D65 white point
def cielab2rgb(L, a, b):
    Ltmp = (L + 16.0) / 116.0
    X = 0.950456 * lab_inverse(Ltmp + a / 500.0)
    Y = 1.000000 * lab_inverse(Ltmp)
    Z = 1.088754 * lab_inverse(Ltmp - b / 200.0)
    Rtmp =  3.2406 * X - 1.5372 * Y - 0.4986 * Z
    Gtmp = -0.9689 * X + 1.8758 * Y + 0.0415 * Z
    Btmp =  0.0557 * X - 0.2040 * Y + 1.0570 * Z
    M = 0.0
    if Rtmp <= Gtmp:
        M = min(Rtmp, Btmp)
    else:
        M = min(Gtmp, Btmp)
    if M < 0:
        Rtmp -= M
        Gtmp -= M
        Btmp -= M
    R = gamma_correction(Rtmp)
    G = gamma_correction(Gtmp)
    B = gamma_correction(Btmp)
    return R, G, B

def bytes2short(a, b, big_endian = False):
    x = 0
    if big_endian == False:
        x = (ord(chr(b)) << 8) + ord(chr(a))
    else:
        x = (ord(chr(a)) << 8) + ord(chr(b))
    return x

def bytes2int(data, big_endian = False):
    ba = bytearray(data)
    if big_endian:
        ba = reversed(ba)
    x = 0
    for offset, byte in enumerate(ba):
        x += byte << (offset * 8)
    return x

def loadSEGmesh(ds):
    ren = vtk.vtkRenderer()
    renwin = vtk.vtkRenderWindow()
    renwin.AddRenderer(ren)
    iren = vtk.vtkRenderWindowInteractor()
    style = vtk.vtkInteractorStyleTrackballCamera()
    iren.SetInteractorStyle(style)
    iren.SetRenderWindow(renwin)
    ren.SetBackground(0.3254, 0.3490, 0.3764)
    renwin.SetSize(800, 600)
    iren.Initialize()
    text1 = vtk.vtkTextActor()
    text1.GetTextProperty().SetFontSize(18)
    text1.GetTextProperty().SetColor(1.0, 1.0, 1.0)
    text1.SetInput(str("Loading ... "))
    text1.SetPosition(100, 100)
    ren.AddActor2D(text1)
    renwin.Render()
    x = 0
    count_points = 0
    count_polys  = 0
    for s in ds.SurfaceSequence:
        mapper = vtk.vtkPolyDataMapper()
        #
        # pydicom 2.x.x
        bytes0 = s.SurfacePointsSequence[0].PointCoordinatesData
        coords = unpack(f"<{int(len(bytes0) / 4)}f", bytes0)
        #
        # pydicom 1.x.x
        #coords = s.SurfacePointsSequence[0].PointCoordinatesData
        #
        points = vtk.vtkPoints()
        for pindx in range(int(len(coords) / 3)):
            point = coords[pindx * 3 : pindx * 3 + 3]
            if dicom_lps_to_slicer_ras:
                point[0] *= -1
                point[1] *= -1
            points.InsertNextPoint(point)
            count_points += 1
        index_type_int = True
        if "LongTrianglePointIndexList" in s.SurfaceMeshPrimitivesSequence[0]:
            tidx = s.SurfaceMeshPrimitivesSequence[0].LongTrianglePointIndexList
        else:
            tidx = s.SurfaceMeshPrimitivesSequence[0].TrianglePointIndexList
            if force_index_type_int == False:
                index_type_int = False
        polys = vtk.vtkCellArray()
        z = 0
        if index_type_int == False:
            while z < len(tidx):
                # 6 bytes to 3 words
                polys.InsertNextCell(3)
                if index_start_with_0 == False:
                    polys.InsertCellPoint(bytes2short(tidx[z],     tidx[z + 1]) - 1)
                    polys.InsertCellPoint(bytes2short(tidx[z + 2], tidx[z + 3]) - 1)
                    polys.InsertCellPoint(bytes2short(tidx[z + 4], tidx[z + 5]) - 1)
                else:
                    polys.InsertCellPoint(bytes2short(tidx[z],     tidx[z + 1]))
                    polys.InsertCellPoint(bytes2short(tidx[z + 2], tidx[z + 3]))
                    polys.InsertCellPoint(bytes2short(tidx[z + 4], tidx[z + 5]))
                count_polys += 1
                z += 6
        else:
            while z < len(tidx):
                # 12 bytes to 3 dwords
                polys.InsertNextCell(3)
                if index_start_with_0 == False:
                    polys.InsertCellPoint(bytes2int(tidx[z     : z +  3]) - 1)
                    polys.InsertCellPoint(bytes2int(tidx[z + 4 : z +  7]) - 1)
                    polys.InsertCellPoint(bytes2int(tidx[z + 8 : z + 11]) - 1)
                else:
                    polys.InsertCellPoint(bytes2int(tidx[z     : z +  3]))
                    polys.InsertCellPoint(bytes2int(tidx[z + 4 : z +  7]))
                    polys.InsertCellPoint(bytes2int(tidx[z + 8 : z + 11]))
                count_polys += 1
                z += 12
        polyData = vtk.vtkPolyData()
        polyData.SetPoints(points)
        polyData.SetPolys(polys)
        gen_normals = vtk.vtkPolyDataNormals()
        gen_normals.SetInputData(polyData)
        gen_normals.ComputePointNormalsOn()
        gen_normals.ComputeCellNormalsOff()
        gen_normals.Update()
        mapper.SetInputData(gen_normals.GetOutput())
        actor = vtk.vtkActor()
        actor.SetMapper(mapper)
        R = 0.8
        G = 0.5
        B = 0.2
        if ("RecommendedDisplayCIELabValue" in s):
            cielab = s.RecommendedDisplayCIELabValue
            L = (cielab[0] / 65535.0) * 100.0
            a = ((cielab[1] - 32896.0) / 65535.0) * 255.0
            b = ((cielab[2] - 32896.0) / 65535.0) * 255.0
            R, G, B = cielab2rgb(L, a, b)
        else:
            if x > 0:
                R = random.randint(0, 255) / 255.0
                G = random.randint(0, 255) / 255.0
                B = random.randint(0, 255) / 255.0
        actor.GetProperty().SetColor(R, G, B)
        ren.AddActor(actor)
        x += 1
        if write_vtk:
            writer_vtk = vtk.vtkPolyDataWriter()
            writer_vtk.SetInputData(polyData)
            writer_vtk.SetFileName(output_prefix + str(x) + '.vtk')
            writer_vtk.Write()
        if write_stl:
            writer_stl.SetFileName(output_prefix + str(x) + '.stl')
            writer_stl.SetInputConnection(polyData)
            writer_stl.SetFileTypeToBinary()
            writer_stl.Write()
    str0 = " surfaces, "
    str1 = " points, "
    str2 = " triangles"
    if x == 1:
        str0 = " surface, "
    if count_points == 1:
        str1 = " point, "
    if count_polys  == 1:
        str2 = " triangle"
    ren.RemoveActor2D(text1)
    text = vtk.vtkTextActor()
    text.GetTextProperty().SetFontSize(16)
    text.GetTextProperty().SetColor(1.0, 1.0, 1.0)
    text.SetInput(str(x) + str0 + str(count_points) + str1 + str(count_polys) + str2)
    text.SetPosition(4, 4)
    ren.AddActor2D(text)
    ren.ResetCamera()
    ren.ResetCameraClippingRange()
    renwin.Render()
    iren.Start()

def load(loadable):
    import pydicom
    print("pydicom ", pydicom.__version__)
    ds = pydicom.read_file(loadable)
    loadSEGmesh(ds)

if __name__ == "__main__":
    if len(sys.argv) == 2:
        if sys.version_info[0] < 3:
            print("Python 3 is required!")
        else:
            filename = os.path.abspath(sys.argv[1])
            load(filename)
    else:
        print("File name is required")