In [2]:
import vtk
import itk
import SimpleITK as sitk
import pydicom as dicom
import nibabel as nib
import matplotlib.pyplot as plt
import numpy as np
import cv2
import skimage
import os
import sys
In [3]:
def read_dicom_series_pydicom(dir_path):
"""
Reads dicom series from directory with pydicom
:param dir_path: string, directory with dicom series
:return: list of dicom images
"""
files = os.listdir(dir_path)
images = []
for f in files:
filepath = os.path.join(dir_path, f)
if os.path.isfile(filepath):
image = dicom.read_file(filepath)
images.append(image)
return images
In [ ]:
def write_dicom_series(images, dir):
"""
Saves dicom series to directory
:param images: list of pydicom images
:param dir: path to directory where it should be saved
:return:
"""
for image in images:
f_name = image.filename.split('/')[-1]
image.save_as(f'{dir}{f_name}')
In [4]:
# read dicoms from path
dicoms = read_dicom_series_pydicom('./3/')
In [5]:
# informations of photo (tags)
print(dicoms[0])
# get numpy array from dicom image
dcm_arr = dicoms[100].pixel_array
# show image in matplotlib
plt.figure(figsize=(10, 10))
plt.imshow(dcm_arr, cmap='gray')
Out[5]:
In [6]:
# read nifti image from directory
filename = os.path.join(nib.testing.data_path, 'example4d.nii.gz')
nii_image = nib.load(filename)
# get numpy ndarray from nifti image
nii_image_arr = nii_image.get_data()
In [7]:
# print information and tags
print(nii_image)
plt.figure(figsize=(10, 10))
plt.imshow(nii_image_arr[:, :, 0, 0], cmap='gray')
Out[7]:
In [8]:
def read_dicom(img_path):
"""
Reads dicom series from directory with SimpleITK
:param img_path: path to directory with a series
:return: list of dicom images
"""
reader = sitk.ImageSeriesReader()
names = reader.GetGDCMSeriesFileNames(img_path)
reader.SetFileNames(names)
reader.MetaDataDictionaryArrayUpdateOn()
reader.LoadPrivateTagsOn()
image = reader.Execute()
return image
In [9]:
# read dicoms and get numpy ndarray from SimpleITK Image
sitk_dicoms = read_dicom('./3/')
sitk_dicom_arr = sitk.GetArrayFromImage(sitk_dicoms)
In [10]:
# print information from SimpleITK Image
print(sitk_dicoms)
plt.figure(figsize=(10, 10))
plt.imshow(sitk_dicom_arr[0], cmap='gray')
Out[10]:
In [5]:
hounsfield_table = cv2.imread('hu_scale_table.png')
plt.figure(figsize=(15, 15))
plt.imshow(hounsfield_table)
def apply_hu(dicoms_list):
"""
Applies Hounsfield units transformation to DICOM images.
Args:
dicoms_list (list): list of DICOM objects (pydicom.dataset.FileDataset).
Returns:
dicoms_hu (list): list of DICOM objects (pydicom.dataset.FileDataset).
"""
dicoms_hu = []
for (i, dcm) in enumerate(dicoms_list):
try:
slope = dcm.RescaleSlope
except Exception as ex:
sys.exit(1)
try:
intercept = dcm.RescaleIntercept
except Exception as ex:
sys.exit(1)
pixels = dicoms_list[i].pixel_array.copy()
pixels = slope * pixels.astype(np.float64)
pixels = pixels.astype(np.int16)
pixels += np.int16(intercept)
dicoms_hu.append(pixels)
return dicoms_hu
In [13]:
# simple operators works on SimpleITK Images as on numpy arrays
s = sitk_dicoms**2
sx = sitk.GetArrayFromImage(s)
x = np.hstack((sx[0], sitk_dicom_arr[0]))
plt.figure(figsize=(20, 10))
plt.imshow(x, cmap='gray')
Out[13]:
In [ ]:
# applies hounsfield scale and threshold (bones and vessels on picture)
hu = apply_hu(dicoms)
threshold = hu.copy()
threshold = np.asarray(threshold)
threshold[threshold > 110] = 255
threshold[threshold <= 110] = 0
plt.figure(figsize=(10, 10))
plt.imshow(threshold[100], cmap='gray')
In [16]:
def show_3d_image_slices(image_array):
"""
Displays 3D image as 2D slices
:key k: left on trackbar
:key l: right on trackbar
:param image_array: 3D image in numpy.ndarray
"""
def nothing(x):
pass
w, h = image_array[0].shape
cv2.namedWindow('name', cv2.WINDOW_NORMAL)
cv2.resizeWindow('name', w, h)
cv2.createTrackbar('image_index', 'name', 0, len(image_array) - 1, nothing)
while 1:
img_index = cv2.getTrackbarPos('image_index', 'name')
show_img = image_array[img_index]
# show_img = show_img.astype(np.uint8)
cv2.imshow('name', show_img)
k = cv2.waitKey(1) & 0xFF
if k == 27:
break
elif k == 107:
if img_index > 0:
cv2.setTrackbarPos('image_index', 'name', img_index - 1)
elif k == 108:
if img_index < len(image_array) - 1:
cv2.setTrackbarPos('image_index', 'name', img_index + 1)
cv2.destroyAllWindows()
In [17]:
# print image in slices
threshold = threshold.astype(np.uint8)
show_3d_image_slices(threshold)
In [16]:
from skimage import measure
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from scipy import ndimage
def plot_3d(image, name):
# p = image.transpose(2, 1, 0)
verts, faces = measure.marching_cubes_classic(image)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
mesh = Poly3DCollection(verts[faces], alpha=0.7)
face_color = [0.45, 0.45, 0.75]
mesh.set_facecolor(face_color)
ax.add_collection(mesh)
ax.set_xlim(0, image.shape[0])
ax.set_ylim(0, image.shape[1])
ax.set_zlim(0, image.shape[2])
plt.savefig(name)
def resample(image, scan, new_spacing=[1, 1, 1]):
# Determine current pixel spacing
spacing = np.array([scan[0].SliceThickness] + [float(x) for x in scan[0].PixelSpacing], dtype=np.float32)
resize_factor = spacing / new_spacing
new_real_shape = image.shape * resize_factor
new_shape = np.round(new_real_shape)
real_resize_factor = new_shape / image.shape
new_spacing = spacing / real_resize_factor
image = ndimage.interpolation.zoom(image, real_resize_factor, mode='nearest')
return image, new_spacing
In [17]:
# plot 3d image in matplotlib (very slow)
resampled, spacing = resample(threshold, dicoms)
plot_3d(resampled, 'matplotlib_3d')
In [8]:
import pclpy
def print_3d_point_cloud(window_name, cloud):
"""
Shows 3D point cloud
:param window_name: string
:param cloud: pclpy.PointCloud.[PointXYZ, PointXYZRGB, ...]
"""
visual = pclpy.pcl.visualization.CloudViewer(window_name)
visual.showCloud(cloud)
v = True
while v:
v = not (visual.wasStopped())
In [9]:
# prepare data to point cloud (works on Windows, no pclpy on Linux)
points = np.nonzero(threshold)
points = np.transpose(points)
point_cloud = pclpy.pcl.PointCloud.PointXYZ(points)
print_3d_point_cloud('', point_cloud)
In [ ]:
import vtk
# vtk example
colors = vtk.vtkNamedColors()
fileName = 'FullHead.mhd' #get_program_parameters()
colors.SetColor("SkinColor", [255, 125, 64, 255])
colors.SetColor("BkgColor", [51, 77, 102, 255])
# Create the renderer, the render window, and the interactor. The renderer
# draws into the render window, the interactor enables mouse- and
# keyboard-based interaction with the data within the render window.
#
aRenderer = vtk.vtkRenderer()
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(aRenderer)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
reader = vtk.vtkMetaImageReader()
reader.SetFileName(fileName)
# An isosurface, or contour value of 500 is known to correspond to the
# skin of the patient.
skinExtractor = vtk.vtkMarchingCubes()
skinExtractor.SetInputConnection(reader.GetOutputPort())
skinExtractor.SetValue(0, 500)
skinMapper = vtk.vtkPolyDataMapper()
skinMapper.SetInputConnection(skinExtractor.GetOutputPort())
skinMapper.ScalarVisibilityOff()
skin = vtk.vtkActor()
skin.SetMapper(skinMapper)
skin.GetProperty().SetDiffuseColor(colors.GetColor3d("SkinColor"))
# An outline provides context around the data.
#
outlineData = vtk.vtkOutlineFilter()
outlineData.SetInputConnection(reader.GetOutputPort())
mapOutline = vtk.vtkPolyDataMapper()
mapOutline.SetInputConnection(outlineData.GetOutputPort())
outline = vtk.vtkActor()
outline.SetMapper(mapOutline)
outline.GetProperty().SetColor(colors.GetColor3d("Black"))
# It is convenient to create an initial view of the data. The FocalPoint
# and Position form a vector direction. Later on (ResetCamera() method)
# this vector is used to position the camera to look at the data in
# this direction.
aCamera = vtk.vtkCamera()
aCamera.SetViewUp(0, 0, -1)
aCamera.SetPosition(0, -1, 0)
aCamera.SetFocalPoint(0, 0, 0)
aCamera.ComputeViewPlaneNormal()
aCamera.Azimuth(30.0)
aCamera.Elevation(30.0)
# Actors are added to the renderer. An initial camera view is created.
# The Dolly() method moves the camera towards the FocalPoint,
# thereby enlarging the image.
aRenderer.AddActor(outline)
aRenderer.AddActor(skin)
aRenderer.SetActiveCamera(aCamera)
aRenderer.ResetCamera()
aCamera.Dolly(1.5)
# Set a background color for the renderer and set the size of the
# render window (expressed in pixels).
aRenderer.SetBackground(colors.GetColor3d("BkgColor"))
renWin.SetSize(640, 480)
# Note that when camera movement occurs (as it does in the Dolly()
# method), the clipping planes often need adjusting. Clipping planes
# consist of two planes: near and far along the view direction. The
# near plane clips out objects in front of the plane the far plane
# clips out objects behind the plane. This way only what is drawn
# between the planes is actually rendered.
aRenderer.ResetCameraClippingRange()
# Initialize the event loop and then start it.
iren.Initialize()
iren.Start()