Skip to content

Instantly share code, notes, and snippets.

@issakomi

issakomi/loadSEGmesh.py

Last active Dec 19, 2020
Embed
What would you like to do?
View DICOM SEG mesh with VTK and pydicom
# 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 bit words,
# but bad datasets exist with 32 bit indices
force_index_type_int = False
#
#
#
#
#
#
#
#
import glob
import os
import sys
import string
import random
import vtk
from struct import unpack
#
#
# CIELab to RGB
# http://www.getreuer.info/home/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):
if isinstance(data, str):
data = bytearray(data)
if big_endian:
data = reversed(data)
x = 0
for offset, byte in enumerate(data):
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")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment