conversion from stl to nifti and vice versa

README.md

• for conversion of nifti to vtk stl, see nifti2stl.py below

• for conversion from stl to nifgi, see stl2nifti.py below

wget https://discourse.itk.org/uploads/default/original/1X/ced4efe8c1eb78342afe474f53ad4409ad5dfc0c.obj
mv ced4efe8c1eb78342afe474f53ad4409ad5dfc0c.obj bunny.obj

wget https://upload.wikimedia.org/wikipedia/commons/4/43/Stanford_Bunny.stl


docker build -t stl2nifti .
docker run -it -u $(id -u):$(id -g) -v $PWD:$PWD -w $PWD stl2nifti bash
python stl2nifti.py bunny.obj bunny.nii.gz

python stl2nifti.py Stanford_Bunny.stl Stanford_Bunny.nii.gz

https://discourse.itk.org/t/trianglemeshtobinaryimagefilter-in-python/1604/7

https://examples.itk.org/src/core/mesh/converttrianglemeshtobinaryimage/documentation

https://github.com/InsightSoftwareConsortium/ITKIOMeshSTL
https://stackoverflow.com/questions/42281881/get-stl-format-3d-mesh-from-binary-mask-segmentation

https://discourse.slicer.org/t/itk-image-to-mesh-to-vtk-polydata/12717

https://slicer.readthedocs.io/en/latest/developer_guide/script_repository.html#rasterize-a-model-and-save-it-to-a-series-of-image-files

https://github.com/InsightSoftwareConsortium/ITK/issues/2416

https://stackoverflow.com/questions/31430171/convert-stl-into-numpy-array

https://github.com/InsightSoftwareConsortium/ITK/issues/2884

https://stackoverflow.com/questions/65112987/convert-stl-object-to-vtk-geometry-in-python/68584066#68584066

DICOM RT-STRUCT

https://zhangresearch.org/post/parse-dicom-rtstruct-into-binary-masks
https://github.com/qurit/rt-utils

Dockerfile

FROM python:3.8-bullseye
COPY requirements.txt /tmp/requirements.txt
RUN pip install -r /tmp/requirements.txt

nifti2stl.py

'''
http://vtk.1045678.n5.nabble.com/How-to-write-vtkDelaunay3D-into-vtk-file-td5741818.html
https://gist.github.com/pangyuteng/facd430d0d9761fc67fff4ff2e5fffc3
'''

import sys
import itk
import numpy as np
import vtk

nifti_file = sys.argv[1]

output_stl_file = sys.argv[2]

reader = vtk.vtkNIFTIImageReader()
reader.SetFileName(nifti_file)
reader.Update()

threshold = vtk.vtkImageThreshold ()
threshold.SetInputConnection(reader.GetOutputPort())
threshold.ThresholdByLower(2000)  #th
threshold.ReplaceInOn()
threshold.SetInValue(0)  # set all values below th to 0
threshold.ReplaceOutOn()
threshold.SetOutValue(1)  # set all values above th to 1
threshold.Update()

dmc = vtk.vtkDiscreteMarchingCubes()
dmc.SetInputConnection(threshold.GetOutputPort())
dmc.GenerateValues(1, 1, 1)
dmc.Update()

mapper = vtk.vtkPolyDataMapper()
mapper.SetInputConnection(dmc.GetOutputPort())
mapper.Update()

smooth = vtk.vtkSmoothPolyDataFilter()
smooth.SetInputConnection(dmc.GetOutputPort())
smooth.SetNumberOfIterations(100)
smooth.SetRelaxationFactor(1)
smooth.Update()

if output_stl_file.endswith('.vtk'):
    writer = vtk.vtkPolyDataWriter()
    writer.SetInputData(smooth.GetOutput())
    writer.SetFileName('bunny.vtk')
    writer.Write()

if output_stl_file.endswith('.ply'):
    writer = vtk.vtkPLYWriter()
    writer.SetInputConnection(dmc.GetOutputPort())
    writer.SetFileName('bunny.ply')
    writer.Write()
    
if output_stl_file.endswith('.stl'):
    writer = vtk.vtkSTLWriter()
    writer.SetInputConnection(dmc.GetOutputPort())
    writer.SetFileTypeToBinary()
    writer.SetFileName("bunny.stl")
    writer.Write()

requirements.txt

itk==5.3.0
itk-iomeshstl==2.3.0

stl2nifti.py

import sys
import itk
import numpy as np

mesh_file = sys.argv[1] # input stl or obj filepath
nifti_file = sys.argv[2] # output nifti filepath

Dimension = 3
PixelType = itk.SS
ImageType = itk.Image[PixelType, Dimension]

meshType = itk.Mesh[itk.SS, 3]
meshReader = itk.MeshFileReader[meshType].New()

if mesh_file.endswith('.obj'):
    meshIO = itk.OBJMeshIO.New()
elif mesh_file.endswith('.stl'):
    meshIO = itk.STLMeshIO.New()
else:
    raise NotImplementedError()
meshReader.SetMeshIO(meshIO)
meshReader.SetFileName(mesh_file)
meshReader.Update()

print(meshReader.GetOutput().GetBoundingBox())

origin = meshReader.GetOutput().GetBoundingBox().GetMinimum()
mesh_size = np.array(meshReader.GetOutput().GetBoundingBox().GetMaximum()) -\
    np.array(meshReader.GetOutput().GetBoundingBox().GetMinimum())

image_size = [512,512,512]
sp = np.max(mesh_size/np.array(image_size))
spacing = [sp]*3
print('mesh_origin',origin)
print('mesh_size',mesh_size)
print('image_size',image_size)
print('spacing',spacing)

region = itk.ImageRegion[Dimension]()
region.SetSize(image_size)
region.SetIndex([0,0,0])
image = ImageType.New()
image.SetRegions(region)
image.Allocate()
image.SetOrigin(origin)
image.SetSpacing(spacing)

myfilter = itk.TriangleMeshToBinaryImageFilter.New()
myfilter.SetInput(meshReader.GetOutput())
myfilter.SetInfoImage(image)
myfilter.SetInsideValue(1)
myfilter.Update()

print(myfilter.GetSize())
print(myfilter.GetSpacing())
print(myfilter.GetOrigin())
print(myfilter.GetDirection())
itk.imwrite(myfilter.GetOutput(),nifti_file)

# https://examples.vtk.org/site/Python/Meshes/Decimation
reduction = 0.9
decimate = vtkDecimatePro()
decimate.SetInputData(dmc.GetOutputPort())
decimate.SetTargetReduction(reduction)
decimate.PreserveTopologyOn()
decimate.Update()
decimated = vtkPolyData()
decimated.ShallowCopy(decimate.GetOutput())

cross references:

nifti to gii
https://github.com/pangyuteng/papaya-flask-celery/blob/master/papaya-gii/nii2gii.py

isosurface rendering with vtk
https://gist.github.com/pangyuteng/facd430d0d9761fc67fff4ff2e5fffc3