Last active
January 13, 2017 05:19
-
-
Save brotherofken/42f3cb84d4e4e5395708f346314240a6 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# In[] | |
import SimpleITK as sitk | |
import vtk | |
import numpy as np | |
import scipy as sp | |
from scipy.stats.mstats import mquantiles | |
import sys | |
from vtk.util.vtkConstants import * | |
##### | |
def numpy2VTK(img,spacing=[1.0,1.0,1.0]): | |
# evolved from code from Stou S., | |
# on http://www.siafoo.net/snippet/314 | |
importer = vtk.vtkImageImport() | |
img_data = img.astype('uint8') | |
img_string = img_data.tostring() # type short | |
dim = img.shape | |
importer.CopyImportVoidPointer(img_string, len(img_string)) | |
importer.SetDataScalarType(VTK_UNSIGNED_CHAR) | |
importer.SetNumberOfScalarComponents(1) | |
extent = importer.GetDataExtent() | |
importer.SetDataExtent(extent[0], extent[0] + dim[2] - 1, | |
extent[2], extent[2] + dim[1] - 1, | |
extent[4], extent[4] + dim[0] - 1) | |
importer.SetWholeExtent(extent[0], extent[0] + dim[2] - 1, | |
extent[2], extent[2] + dim[1] - 1, | |
extent[4], extent[4] + dim[0] - 1) | |
importer.SetDataSpacing( spacing[0], spacing[1], spacing[2]) | |
importer.SetDataOrigin( 0,0,0 ) | |
return importer | |
def volumeRender(img, tf=[],spacing=[1.0,1.0,1.0]): | |
importer = numpy2VTK(img,spacing) | |
# Transfer Functions | |
opacity_tf = vtk.vtkPiecewiseFunction() | |
color_tf = vtk.vtkColorTransferFunction() | |
if len(tf) == 0: | |
tf.append([img.min(),0,0,0,0]) | |
tf.append([img.max(),1,1,1,1]) | |
for p in tf: | |
color_tf.AddRGBPoint(p[0], p[1], p[2], p[3]) | |
opacity_tf.AddPoint(p[0], p[4]) | |
# working on the GPU | |
volMapper = vtk.vtkGPUVolumeRayCastMapper() | |
volMapper.SetInputConnection(importer.GetOutputPort()) | |
# The property describes how the data will look | |
volProperty = vtk.vtkVolumeProperty() | |
volProperty.SetColor(color_tf) | |
volProperty.SetScalarOpacity(opacity_tf) | |
volProperty.ShadeOn() | |
volProperty.SetInterpolationTypeToLinear() | |
# # working on the CPU | |
# volMapper = vtk.vtkVolumeRayCastMapper() | |
# compositeFunction = vtk.vtkVolumeRayCastCompositeFunction() | |
# compositeFunction.SetCompositeMethodToInterpolateFirst() | |
# volMapper.SetVolumeRayCastFunction(compositeFunction) | |
# volMapper.SetInputConnection(importer.GetOutputPort()) | |
# | |
# # The property describes how the data will look | |
# volProperty = vtk.vtkVolumeProperty() | |
# volProperty.SetColor(color_tf) | |
# volProperty.SetScalarOpacity(opacity_tf) | |
# volProperty.ShadeOn() | |
# volProperty.SetInterpolationTypeToLinear() | |
# Do the lines below speed things up? | |
# pix_diag = 5.0 | |
# volMapper.SetSampleDistance(pix_diag / 5.0) | |
# volProperty.SetScalarOpacityUnitDistance(pix_diag) | |
vol = vtk.vtkVolume() | |
vol.SetMapper(volMapper) | |
vol.SetProperty(volProperty) | |
return [vol] | |
def vtk_basic( actors ): | |
""" | |
Create a window, renderer, interactor, add the actors and start the thing | |
Parameters | |
---------- | |
actors : list of vtkActors | |
Returns | |
------- | |
nothing | |
""" | |
# create a rendering window and renderer | |
ren = vtk.vtkRenderer() | |
renWin = vtk.vtkRenderWindow() | |
renWin.AddRenderer(ren) | |
renWin.SetSize(600,600) | |
# ren.SetBackground( 1, 1, 1) | |
# create a renderwindowinteractor | |
iren = vtk.vtkRenderWindowInteractor() | |
iren.SetRenderWindow(renWin) | |
for a in actors: | |
# assign actor to the renderer | |
ren.AddActor(a ) | |
# render | |
renWin.Render() | |
# enable user interface interactor | |
iren.Initialize() | |
iren.Start() | |
# In[] Load series | |
from vtk.util import numpy_support | |
path = "/home/rakhunzy/lung_samples/0a0c32c9e08cc2ea76a71649de56be6d" | |
PathDicom = path | |
reader = vtk.vtkDICOMImageReader() | |
reader.SetDirectoryName(PathDicom) | |
reader.Update() | |
_extent = reader.GetDataExtent() | |
ConstPixelDims = [_extent[1]-_extent[0]+1, _extent[3]-_extent[2]+1, _extent[5]-_extent[4]+1] | |
ConstPixelSpacing = reader.GetPixelSpacing() | |
def vtkImageToNumPy(image, pixelDims): | |
pointData = image.GetPointData() | |
arrayData = pointData.GetArray(0) | |
ArrayDicom = numpy_support.vtk_to_numpy(arrayData) | |
ArrayDicom = ArrayDicom.reshape(pixelDims, order='F') | |
return ArrayDicom | |
reader_data = vtkImageToNumPy(reader.GetOutput(), ConstPixelDims).astype(np.float) | |
image_data = 255 * (reader_data - reader_data.min()) / (reader_data.max() - reader_data.min()) | |
_extent = reader.GetDataExtent() | |
ConstPixelDims = [_extent[1]-_extent[0]+1, _extent[3]-_extent[2]+1, _extent[5]-_extent[4]+1] | |
ConstPixelSpacing = reader.GetPixelSpacing() | |
# In[] Threshold to see bones | |
threshold = vtk.vtkImageThreshold () | |
threshold.SetInputConnection(reader.GetOutputPort()) | |
threshold.ThresholdByLower(400) # remove all soft tissue | |
threshold.ReplaceInOn() | |
threshold.SetInValue(0) # set all values below 400 to 0 | |
threshold.ReplaceOutOn() | |
threshold.SetOutValue(1) # set all values above 400 to 1 | |
threshold.Update() | |
# In[] Rescale data | |
threshold_data = vtkImageToNumPy(threshold.GetOutput(), ConstPixelDims).astype(np.float) | |
image_data = 255 * (threshold_data - threshold_data.min()) / (threshold_data.max() - threshold_data.min()) | |
# In[] Render in window | |
from scipy.stats.mstats import mquantiles | |
q = mquantiles(image_data.flatten(),[0.7,0.98]) | |
q[0]=max(q[0],1) | |
q[1] = max(q[1],1) | |
tf=[[0,0,0,0,0],[q[0],0,0,0,0],[q[1],1,1,1,0.5],[image_data.max(),1,1,1,1]] | |
actor_list = volumeRender(image_data, tf=tf, spacing=ConstPixelSpacing) | |
vtk_basic(actor_list) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment