import os
import numpy
import SimpleITK
import matplotlib.pyplot as plt
%pylab inline
Populating the interactive namespace from numpy and matplotlib
We're gonna use this function to quickly plot a 2D SimpleITK image with a grayscale colormap and accompanying axes
def sitk_show(img, title=None, margin=0.05, dpi=40 ):
nda = SimpleITK.GetArrayFromImage(img)
spacing = img.GetSpacing()
figsize = (1 + margin) * nda.shape[0] / dpi, (1 + margin) * nda.shape[1] / dpi
extent = (0, nda.shape[1]*spacing[1], nda.shape[0]*spacing[0], 0)
fig = plt.figure(figsize=figsize, dpi=dpi)
ax = fig.add_axes([margin, margin, 1 - 2*margin, 1 - 2*margin])
plt.set_cmap("gray")
ax.imshow(nda,extent=extent,interpolation=None)
if title:
plt.title(title)
plt.show()
# Directory where the DICOM files are being stored (in this
# case the 'MyHead' folder).
pathDicom = "./MyHead/"
# Z slice of the DICOM files to process. In the interest of
# simplicity, segmentation will be limited to a single 2D
# image but all processes are entirely applicable to the 3D image
idxSlice = 50
# int labels to assign to the segmented white and gray matter.
# These need to be different integers but their values themselves
# don't matter
labelWhiteMatter = 1
labelGrayMatter = 2
Load and read-in the DICOM files
reader = SimpleITK.ImageSeriesReader()
filenamesDICOM = reader.GetGDCMSeriesFileNames(pathDicom)
reader.SetFileNames(filenamesDICOM)
imgOriginal = reader.Execute()
Limit image to a single 2D slice
imgOriginal = imgOriginal[:,:,idxSlice]
Visualize original data
sitk_show(imgOriginal)
Apply a curvature-flow filter to smoothen the image (get rid of image noise)
imgSmooth = SimpleITK.CurvatureFlow(image1=imgOriginal,
timeStep=0.125,
numberOfIterations=5)
# blurFilter = SimpleITK.CurvatureFlowImageFilter()
# blurFilter.SetNumberOfIterations(5)
# blurFilter.SetTimeStep(0.125)
# imgSmooth = blurFilter.Execute(imgOriginal)
sitk_show(imgSmooth)
ConnectedThreshold
filter¶Initial segmentation of white matter
lstSeeds = [(150,75)]
imgWhiteMatter = SimpleITK.ConnectedThreshold(image1=imgSmooth,
seedList=lstSeeds,
lower=130,
upper=190,
replaceValue=labelWhiteMatter)
Visualization with label overlay
# Rescale 'imgSmooth' and cast it to an integer type to match that of 'imgWhiteMatter'
imgSmoothInt = SimpleITK.Cast(SimpleITK.RescaleIntensity(imgSmooth), imgWhiteMatter.GetPixelID())
# Use 'LabelOverlay' to overlay 'imgSmooth' and 'imgWhiteMatter'
sitk_show(SimpleITK.LabelOverlay(imgSmoothInt, imgWhiteMatter))
Fill holes in the white matter label through VotingBinaryHoleFilling
imgWhiteMatterNoHoles = SimpleITK.VotingBinaryHoleFilling(image1=imgWhiteMatter,
radius=[2]*3,
majorityThreshold=1,
backgroundValue=0,
foregroundValue=labelWhiteMatter)
sitk_show(SimpleITK.LabelOverlay(imgSmoothInt, imgWhiteMatterNoHoles))
Segmentation & hole-filling of gray matter
lstSeeds = [(119, 83), (198, 80), (185, 102), (164, 43)]
imgGrayMatter = SimpleITK.ConnectedThreshold(image1=imgSmooth,
seedList=lstSeeds,
lower=150,
upper=270,
replaceValue=labelGrayMatter)
imgGrayMatterNoHoles = SimpleITK.VotingBinaryHoleFilling(image1=imgGrayMatter,
radius=[2]*3,
majorityThreshold=1,
backgroundValue=0,
foregroundValue=labelGrayMatter)
sitk_show(SimpleITK.LabelOverlay(imgSmoothInt, imgGrayMatterNoHoles))
Bitwire OR combination of the white matter and gray matter labels
imgLabels = imgWhiteMatterNoHoles | imgGrayMatterNoHoles
sitk_show(SimpleITK.LabelOverlay(imgSmoothInt, imgLabels))
Assignment of common areas to the gray matter
imgGrayMatterNoHoles.
2
SimpleITK.GetArrayFromImage(imgWhiteMatterNoHoles).min()
0
# imgMask= imgWhiteMatterNoHoles/labelWhiteMatter * imgGrayMatterNoHoles/labelGrayMatter
imgMask= imgWhiteMatterNoHoles * imgGrayMatterNoHoles
imgWhiteMatterNoHoles -= 1000*imgMask*labelWhiteMatter
imgLabels = imgWhiteMatterNoHoles + imgGrayMatterNoHoles
sitk_show(SimpleITK.LabelOverlay(imgSmoothInt, imgLabels))
sitk_show(SimpleITK.LabelOverlay(imgSmoothInt, SimpleITK.LabelContour(imgLabels)))