Python version: 3.4.3 (Anaconda 2.3.0)
IPython version: 4.0.0
SimpleITK version: 0.9.0-g45317
Operating system: OS X El Capitan v10.11.1
This medical imaging project is hosted here on Github.
In [1]:
import collections
In [2]:
import matplotlib
In [3]:
matplotlib.use('Agg')
In [26]:
%matplotlib inline
In [5]:
import importlib
In [6]:
import matplotlib.pyplot as plt
In [7]:
import numpy as np
In [8]:
import scipy.stats
In [9]:
import SimpleITK as sitk
In [10]:
from os.path import expanduser, join
In [11]:
from scipy.spatial.distance import euclidean
In [12]:
import imgscroll
In [13]:
def show_imgs(imgstack, imgstack2=None):
"""Uses imgscroll to show stack of images in a new window with
mouse scrolling. If a 2nd image stack is provided, then it will
be overlaid on top of the 1st. %matplotlib qt must used!"""
X = sitk.GetArrayFromImage(imgstack)
fig = plt.figure()
ax = fig.add_subplot(111)
if imgstack2 is not None:
X2 = sitk.GetArrayFromImage(imgstack2)
maskVal = scipy.stats.mode(X2.flatten())[0][0]
Xmask = np.ma.masked_where(X2 == maskVal, X2)
tracker = imgscroll.IndexTracker2(ax, X, Xmask)
else:
tracker = imgscroll.IndexTracker(ax, X)
fig.canvas.mpl_connect('scroll_event', tracker.onscroll)
plt.show()
In [14]:
def input_level_set_click(featImg, slice2coords):
"""Create input level set from coordinates gathered from user mouse clicks.
A circle of size RADIUS is drawn around each coordinate."""
RADIUS = 10
numCols = featImg.GetSize()[0]
numRows = featImg.GetSize()[1]
depthSize = featImg.GetSize()[2]
X = np.zeros((depthSize, numRows, numCols), dtype=np.int)
# draw circle by setting points within radius of coordinate equal to 1
for n in slice2coords:
for c in slice2coords[n]:
rowIni, rowEnd = c[1] - RADIUS, c[1] + RADIUS
colIni, colEnd = c[0] - RADIUS, c[0] + RADIUS
for i in range(rowIni, rowEnd+1):
for j in range(colIni, colEnd+1):
if euclidean((i,j), (c[1], c[0])) <= RADIUS:
X[n,i,j] = 1
dilation = sitk.BinaryDilate(sitk.GetImageFromArray(X), 3)
img = sitk.Cast(dilation, featImg.GetPixelIDValue()) * -1 + 0.5
img.CopyInformation(featImg)
return img
In [15]:
class IndexMouseCapture(object):
"""For a stack of images, capture and save coordinates from mouse clicks.
A circle of size RADIUS will be displayed upon each click."""
def __init__(self, ax, X):
self.ax = ax
self.X = X
self.slices, numrows, numcols = X.shape
self.ind = 0
self.ind2coord = collections.defaultdict(list)
self.circles = list()
self.im = ax.imshow(self.X[self.ind, :, :], cmap=plt.cm.Greys_r)
self.update()
def onscroll(self, event):
if event.button == 'up':
self.ind = np.clip(self.ind+1, 0, self.slices-1)
else:
self.ind = np.clip(self.ind-1, 0, self.slices-1)
for circ in self.circles:
circ.remove()
self.circles = list()
self.update()
def onclick(self, event):
RADIUS = 10
ix, iy = int(round(event.xdata)), int(round(event.ydata))
self.ind2coord[self.ind].append((ix, iy))
circ = plt.Circle((ix, iy), RADIUS, color='b')
plt.gcf().gca().add_artist(circ)
self.circles.append(circ)
self.im.axes.figure.canvas.draw()
def update(self):
self.im.set_data(self.X[self.ind, :, :])
self.ax.set_title('slice %s' %self.ind)
self.im.axes.figure.canvas.draw()
In [16]:
def tile_slices(imgstack, slicesToShow, imgstack2=None):
"""Show a number of slices tiled side-by-side"""
fig = plt.figure(figsize=(15,10))
numSlices = len(slicesToShow)
for i,sliceNum in enumerate(slicesToShow):
plt.subplot(1, numSlices, i+1)
plt.imshow(sitk.GetArrayFromImage(imgstack[:,:,sliceNum]), cmap=plt.cm.Greys_r)
if imgstack2 is not None:
X2 = sitk.GetArrayFromImage(imgstack2[:,:,sliceNum])
maskVal = scipy.stats.mode(X2.flatten())[0][0]
Xmask = np.ma.masked_where(X2 == maskVal, X2)
plt.imshow(Xmask, alpha=0.5)
plt.gca().get_xaxis().set_ticks([])
plt.gca().get_yaxis().set_ticks([])
fig.subplots_adjust(hspace=0.0, wspace=0.0)
plt.show()
The case being analyzed here is a resectable hepatocellular carcinoma from The Cancer Imaging Archive. The patient ID is TCGA-BC-4073. Ultimately, the liver along with the tumor and blood vessels are to be reconstructed. Here, the focus is only on the liver.
In [17]:
dicomPath = join(expanduser('~'), 'Documents', 'SlicerDICOMDatabase', 'TCIALocal', '0', 'images', '')
reader = sitk.ImageSeriesReader()
seriesIDread = reader.GetGDCMSeriesIDs(dicomPath)[1]
dicomFilenames = reader.GetGDCMSeriesFileNames(dicomPath, seriesIDread)
reader.SetFileNames(dicomFilenames)
imgSeries = reader.Execute()
In [18]:
# choose a series of slices
imgSlices = imgSeries[:,:,30:40]
In [25]:
# %matplotlib qt
# displays scrollable image series in a new window
show_imgs(imgSlices)
Show a few slices, for example, 30, 35, and 39:
In [158]:
tile_slices(imgSlices, [0, 5, 9])
In [19]:
timeStep_, conduct, numIter = (0.04, 9.0, 5)
imgRecast = sitk.Cast(imgSlices, sitk.sitkFloat32)
curvDiff = sitk.CurvatureAnisotropicDiffusionImageFilter()
curvDiff.SetTimeStep(timeStep_)
curvDiff.SetConductanceParameter(conduct)
curvDiff.SetNumberOfIterations(numIter)
imgFilter = curvDiff.Execute(imgRecast)
In [20]:
sigma_ = 2.0
imgGauss = sitk.GradientMagnitudeRecursiveGaussian(image1=imgFilter, sigma=sigma_)
In [20]:
# %matplotlib qt
# display in new window to find values on boundary
show_imgs(imgGauss)
Following Section 4.3.1 "Fast Marching Segmentation" on pages 373-374 from The ITK Software Guide Book 2: Design and Functionality, 4th ed. for the setup of the sigmoid filter. The output will be passed along to a segmentation algorithm below.
Note that the plan is to conduct multiple rounds of segmentation, to "start with a downsampled volume and work back to the full resolution using the results at each intermediate scale as the initialization for the next scale." (pg. 370) Therefore, K1 and K2 below for the sigmoid mapping are first set loosely set and will become more strict in later segmentations.
In [21]:
K1, K2 = 18.0, 8.0
In [22]:
alpha_ = (K2 - K1)/6
beta_ = (K1 + K2)/2
sigFilt = sitk.SigmoidImageFilter()
sigFilt.SetAlpha(alpha_)
sigFilt.SetBeta(beta_)
sigFilt.SetOutputMaximum(1.0)
sigFilt.SetOutputMinimum(0.0)
imgSigmoid = sigFilt.Execute(imgGauss)
In [35]:
# %matplotlib qt
# examine image series in a new window
show_imgs(imgSigmoid)
Show slices 30, 35, and 39 having computed the gradient and sigmoid mapping:
In [26]:
tile_slices(imgSigmoid, [0, 5, 9])
Using ideas from the SimpleITK geodesic active contour example to create an initial input level set. Instead of computing a signed Maurer distance map and then applying a binary threshold, the approach here simply draws a circle of a given radius around each user-chosen seed coordinate. Following the SimpleITK Notebook on Levelset Segmentation, a binary dilation with a kernel size of 3 is performed. Finally, as was done in the example (line 60), all image values are multiplied by -1 and added to 0.5. The results is the input level set.
Use the class IndexMouseCapture (defined above) to capture coordinates from mouse clicks for seeds. The radii are all assumed to be of the same size at the moment.
In [24]:
%matplotlib qt
# get coordinates from mouse clicks (opens in new window)
X = sitk.GetArrayFromImage(imgSigmoid)
fig = plt.figure()
ax = fig.add_subplot(111)
tracker = IndexMouseCapture(ax, X)
fig.canvas.mpl_connect('scroll_event', tracker.onscroll)
fig.canvas.mpl_connect('button_press_event', tracker.onclick)
Out[24]:
In [25]:
# call function to create input level set from the coordinates chosen above by user
initImg = input_level_set_click(imgSigmoid, tracker.ind2coord)
Display slices with the input level set overlaying the sigmoid-mapped image:
In [27]:
tile_slices(imgSigmoid, [4, 5, 6], initImg)
In [28]:
tile_slices(imgSigmoid, [7, 8, 9], initImg)
In [29]:
gac = sitk.GeodesicActiveContourLevelSetImageFilter()
gac.SetPropagationScaling(1.0)
gac.SetCurvatureScaling(0.2)
gac.SetAdvectionScaling(3.0)
gac.SetMaximumRMSError(0.01)
gac.SetNumberOfIterations(200)
Out[29]:
In [30]:
gac3D = gac.Execute(initImg, sitk.Cast(imgSlices, sitk.sitkFloat32))
Displaying the segmentation result from some slices:
In [31]:
tile_slices(imgSigmoid, [4, 5, 6], gac3D)
In [32]:
tile_slices(imgSigmoid, [7, 8, 9], gac3D)
In [ ]: