render isosurface and save as gif using vtk and moviepy

gif_isosurface.md

source /anaconda3/bin/activate py36

wget https://github.com/jonclayden/RNifti/raw/master/inst/extdata/example.nii.gz

python isosurface_vtk.py

isosurface_vtk.py

#
import sys
import numpy as np
import vtk
from vtk.util.colors import salmon


reader = vtk.vtkNIFTIImageReader()
reader.SetFileName('example.nii.gz')
reader.Update()

threshold = vtk.vtkImageThreshold ()
threshold.SetInputConnection(reader.GetOutputPort())
#threshold.ThresholdBetween(1,1)  #th
threshold.ThresholdByLower(50)  #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()

'''
voi = vtk.vtkExtractVOI()
voi.SetInputConnection(threshold.GetOutputPort()) 
voi.SetVOI(0,95, 0,95, 0,59)
voi.SetSampleRate(1,1,1)
#voi.SetSampleRate(3,3,3)
voi.Update()#necessary for GetScalarRange()
srange= voi.GetOutput().GetScalarRange()#needs Update() before!
print("Range", srange)
'''

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

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

actor = vtk.vtkActor()
actor.SetMapper(mapper)
mapper.ScalarVisibilityOff()

actor.GetProperty().SetColor(salmon)
actor.GetProperty().SetOpacity(0.5)
actor.RotateX(90)

renderer = vtk.vtkRenderer()
renderWindow = vtk.vtkRenderWindow()
renderWindow.SetOffScreenRendering(1)
renderWindow.AddRenderer(renderer)

renderer.AddActor(actor)
renderer.SetBackground(1.0, 1.0, 1.0)

center = actor.GetCenter()

camera = renderer.MakeCamera()
camera.SetPosition(0,0,800)
camera.SetFocalPoint(center)

renderer.SetActiveCamera(camera)
renderWindow.Render()

focal_point = camera.GetFocalPoint()
view_up = camera.GetViewUp()
position = camera.GetPosition() 

axis = [0,0,0]
axis[0] = -1*camera.GetViewTransformMatrix().GetElement(0,0)
axis[1] = -1*camera.GetViewTransformMatrix().GetElement(0,1)
axis[2] = -1*camera.GetViewTransformMatrix().GetElement(0,2)

print(position,focal_point,view_up,)

print(camera.GetViewTransformMatrix())
print(camera.GetViewTransformMatrix().GetElement(0,0))
print(camera.GetViewTransformMatrix().GetElement(0,1))
print(camera.GetViewTransformMatrix().GetElement(0,2))

frame_list = []
for n,q in enumerate([10]*35):

    transform = vtk.vtkTransform()
    transform.Identity()

    transform.Translate(*center)
    transform.RotateWXYZ(q,view_up)
    transform.RotateWXYZ(0,axis)
    transform.Translate(*[-1*x for x in center])

    new_position = [0,0,0]
    new_focal_point = [0,0,0]
    transform.TransformPoint(position,new_position)
    transform.TransformPoint(focal_point,new_focal_point)

    camera.SetPosition(new_position)
    camera.SetFocalPoint(new_focal_point)

    focal_point = camera.GetFocalPoint()
    view_up = camera.GetViewUp()
    position = camera.GetPosition() 

    camera.OrthogonalizeViewUp();
    renderer.ResetCameraClippingRange();
    
    renderWindow.Render()
    windowToImageFilter = vtk.vtkWindowToImageFilter()
    windowToImageFilter.SetInput(renderWindow)
    windowToImageFilter.Update()

    writer = vtk.vtkPNGWriter()
    
    fpath = "zisosurface{}.png".format(n)
    writer.SetFileName(fpath)
    writer.SetInputConnection(windowToImageFilter.GetOutputPort())
    writer.Write()

    frame_list.append(fpath)
    
duration = 2
fps = 17.5
time_list = list(np.arange(0,duration,1./fps))
img_dict = {a:f for a,f in zip(time_list,frame_list)}

import PIL
from moviepy import editor
def make_frame(t):    
    fpath= img_dict[t]
    im = PIL.Image.open(fpath)
    ar = np.asarray(im)
    return ar
    
gif_path = 'ani.gif'
clip = editor.VideoClip(make_frame, duration=duration)
clip.write_gif(gif_path, fps=fps)

    
'''
ref
https://www.programcreek.com/python/example/9480/vtk.vtkRenderWindows
https://pyscience.wordpress.com/2014/09/11/surface-extraction-creating-a-mesh-from-pixel-data-using-python-and-vtk/

https://stackoverflow.com/questions/43046798/how-to-rotate-a-vtk-scene-around-a-specific-point
https://gamedev.stackexchange.com/questions/20758/how-can-i-orbit-a-camera-about-its-target-point
'''

save from nifti to stl or vtk files
https://gist.github.com/pangyuteng/0493816e0398cdb59fde48a946745b31
nifti to gii
https://github.com/pangyuteng/papaya-flask-celery/blob/master/papaya-gii/nii2gii.py

# https://gist.github.com/gattia/dc69d6aed71bb6671259
# https://discourse.vtk.org/t/vtk-has-no-attribute-vtk-volume-ray-cast-composite-function/12870/2

# Transfer Functions
opacity_tf = vtk.vtkPiecewiseFunction()
color_tf = vtk.vtkColorTransferFunction()

tf = []
tf.append([0.0,0,0,0,0.0])
tf.append([0.25,0,0,0,0.0])
tf.append([0.25,0,0,0,0.1])
tf.append([1.5,1,0,0,0.2])
tf.append([3.5,1,1,0,0.3])
tf.append([5.0,1,1,1,0.5])

for p in tf:
  color_tf.AddRGBPoint(p[0], p[1], p[2], p[3])
  opacity_tf.AddPoint(p[0], p[4])

volMapper = vtk.vtkOpenGLGPUVolumeRayCastMapper()
volMapper.SetBlendModeToComposite()
volMapper.SetInputConnection(reader.GetOutputPort())

# The property describes how the data will look
volProperty =  vtk.vtkVolumeProperty()
volProperty.SetColor(color_tf)
volProperty.SetScalarOpacity(opacity_tf)
volProperty.ShadeOn()
volProperty.SetInterpolationTypeToLinear()

vol = vtk.vtkVolume()
vol.SetMapper(volMapper)
vol.SetProperty(volProperty)
vol.RotateX(90)

mylist.append(vol)

renderer.AddActor(vol)

vtk-osmesa (no building required!)

https://github.com/pangyuteng/public-scratch/blob/master/docker/dcm/Dockerfile