Segment_all_patients

This is a script to segment all the patients using the procedure outlined in the previous script

BMI 260 Assignment 1

Author: Yusuf Roohani (yroohani@stanford)



In [1]:

    
import dicom
import numpy as np
import cv2

import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
%matplotlib inline

from skimage.filters import threshold_otsu
from skimage.measure import label, regionprops
from skimage.color import label2rgb

from os import listdir
from os.path import join

from mpl_toolkits.mplot3d import Axes3D



In [2]:

    
# This is the functionalized segmentation procedure post-thresholding
# Please add to your path
from Segmenter import vol_segment



In [3]:

    
# List all patients
patients = listdir('./kaggle_sample_data/')

Perform segmentation



In [11]:

    
segmented_lung = {}
vols = {}

# Cycle over all patients
# Takes a while...
for patient in patients:
    
    if patient[0] != '0' or patient == '0b20184e0cd497028bdd155d9fb42dc9':
        continue

    # Read in filenames
    path_dcms = './kaggle_sample_data/' + patient + '/'
    list_dcms = listdir(path_dcms)

    # Create variables to hold a 3D volume
    width = 512; height = 512; depth = len(list_dcms)
    lung_vol_houns = np.zeros([width,height,depth])

    # Create a dictionary to link slice location with image
    dcm_slices_houns = {}

    # Stack up all the slices in order
    for d in range(depth):
        dcm_file = dicom.read_file(path_dcms + list_dcms[d])

        # This is to club bone and other organs into a single broad non-lung class
        dcm_file.pixel_array_houns = (dcm_file.RescaleSlope * dcm_file.pixel_array) + dcm_file.RescaleIntercept
        dcm_file.pixel_array_houns[dcm_file.pixel_array_houns > -450] = -450
        dcm_slices_houns[dcm_file.SliceLocation] = dcm_file.pixel_array_houns

    sorted_slices = np.sort(dcm_slices_houns.keys())

    # Now let's make a single volume
    for idx,s in enumerate(sorted_slices):
        lung_vol_houns[:,:,idx] = dcm_slices_houns[s]

    # Let's max min normalize these images
    lung_vol_norm = (lung_vol_houns - np.min(lung_vol_houns))/(np.max(lung_vol_houns) - np.min(lung_vol_houns))

    # Otsu thresholding
    flat_volume = lung_vol_norm.flatten()
    noise_thresh = 0.05
    flat_volume_clean = flat_volume[flat_volume > noise_thresh]
    otsu_thresh = threshold_otsu(flat_volume_clean)

    # Do the actual segmentation
    segmented_lung[patient], vols[patient]  = vol_segment(lung_vol_norm, otsu_thresh)

Study region volumes to design threshold



In [21]:

    
# Keeping track of the region volumes
bads = []; goods = []
for vol in vols.itervalues():
    bads.append(vol[2])
    goods.append(vol[1])



In [29]:

    
plt.plot(range(len(bads)),bads)
plt.plot(range(len(goods)),goods)
plt.plot([0, len(goods)], [1e6, 1e6], 'k-', lw=2)
ax.lines









    Out[29]:





[]

Create 3D visualizations for the report



In [12]:

    
# Create a 3D visualization
    
for patient in segmented_lung.iterkeys():
    segmented = segmented_lung[patient]

    fig = plt.figure(figsize = (15,15))
    ax = fig.gca(projection='3d')
    ax.set_zlim([0,segmented.shape[2]+4])
    X = np.tile(np.arange(512),[512,1])
    Y = np.transpose(np.tile(np.arange(512),[512,1]))

    for i in range(segmented.shape[2]):
        ax.contour(X, Y, segmented[:,:,i]*(i+1), levels = [ i])

    fig.savefig('./3D1_' + patient + '.png')