Tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# U-Net Segmentation Approach to Cancer Diagnosis\n",
    "*written by Jonathan Mulholland, Data Scientist, Booz Allen Hamilton*\n",
    "\n",
    "In this tutorial, we will show one approach to predicting whether a CT scan is of a patient who either has or will develop cancer within the next 12 months or not.\n",
    "\n",
    "We will train a network to segment out potentially cancerous nodules and then use the characteristics of that segmentation to make predictions about the diagnosis of the scanned patient within a 12 month time frame. \n",
    "\n",
    "The code for this tutorial can be found at https://github.com/booz-allen-hamilton/DSB3Tutorial \n",
    "## Dependencies and tools\n",
    "This tutorial uses python and has the following dependencies\n",
    "\n",
    "* numpy\n",
    "* scikit-image\n",
    "* scikit-learn \n",
    "* keras (tensorflow backend)\n",
    "* matplotlib\n",
    "* pydicom\n",
    "* SimpleITK\n",
    "\n",
    "Note: Keras allows multiple backends for use when training. We chose a GPU enabled installation of tensorflow as the backend of Keras.\n",
    "\n",
    "In order to identify regions with nodules, we will use a U-Net style convolutional network which was designed for segmenting neuronal structures. https://arxiv.org/abs/1505.04597 \n",
    "\n",
    "Our code for the network was based on a turorial posted by Marko Jocic on the Kaggle forum for the Ultrasound Nerve Segmentation Challenge. https://www.kaggle.com/c/ultrasound-nerve-segmentation/forums/t/21358/0-57-deep-learning-keras-tutorial\n",
    "\n",
    "The images we'll be predicting cancer diagnoses on are scans from low-dose helical computed tomography (CT). The appearance on nodules within the CT scan indicate the possibility of cancer, and we will need training examples with marked nodules in order train the U Net to find these nodules. Rather than hand label images, we turn to the Lung Nodule Analysis 2016 (LUNA2016) challenge which has made available CT images with annotated nodule locations. We will first use the LUNA data set to generate an appropriate training set for our U-Net. We will use these examples to train our supervised segmenter. \n",
    "\n",
    "## Constructing a training set from the LUNA 2016 data\n",
    "\n",
    "\n",
    "We are going to use the nodule locations as given in annotations.csv and extract three transverse slices that contain the largest nodule from each patient scan. Masks will be created for those slices based on the nodule dimensions given in annotations.csv. The output of this file will be two files for each patient scan: a set of images and a set of corresponding nodule masks. The data from the LUNA 2016 challenge can be found at https://luna16.grand-challenge.org/\n",
    "\n",
    "First we import the necessary tools and find the largest nodule in the patient scan. There are multiple nodule listings for some patients in annotations.csv. We're using a pandas DataFrame named ```df_node``` to keep track of the case numbers and the node information. The node information is an (x,y,z) coordinate in mm using a coordinate system defined in the .mhd file. \n",
    "\n",
    "The following snippets of code are from LUNA_mask_extraction.py:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "import SimpleITK as sitk\n",
    "import numpy as np\n",
    "import csv\n",
    "from glob import glob\n",
    "import pandas as pd\n",
    "\n",
    "file_list=glob(luna_subset_path+\"*.mhd\")\n",
    "\n",
    "#####################\n",
    "#\n",
    "# Helper function to get rows in data frame associated \n",
    "# with each file\n",
    "def get_filename(case):\n",
    "    global file_list\n",
    "    for f in file_list:\n",
    "        if case in f:\n",
    "            return(f)\n",
    "#\n",
    "# The locations of the nodes\n",
    "df_node = pd.read_csv(luna_path+\"annotations.csv\")\n",
    "df_node[\"file\"] = df_node[\"seriesuid\"].apply(get_filename)\n",
    "df_node = df_node.dropna()\n",
    "\n",
    "#####\n",
    "#\n",
    "# Looping over the image files\n",
    "#\n",
    "fcount = 0\n",
    "for img_file in file_list:\n",
    "    print \"Getting mask for image file %s\" % img_file.replace(luna_subset_path,\"\")\n",
    "    mini_df = df_node[df_node[\"file\"]==img_file] #get all nodules associate with file\n",
    "    if len(mini_df)>0:       # some files may not have a nodule--skipping those \n",
    "        biggest_node = np.argsort(mini_df[\"diameter_mm\"].values)[-1]   # just using the biggest node\n",
    "        node_x = mini_df[\"coordX\"].values[biggest_node]\n",
    "        node_y = mini_df[\"coordY\"].values[biggest_node]\n",
    "        node_z = mini_df[\"coordZ\"].values[biggest_node]\n",
    "        diam = mini_df[\"diameter_mm\"].values[biggest_node]\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting the nodule position in the mhd files\n",
    "\n",
    "The nodule locations are given in terms of millimeters relative to a coordinate system defined by the CT scanner. The image data is given as a varying length stack of 512 X 512 arrays. In order to translate the voxel position to the world coordinate system, one needs to know the real world position of the [0,0,0] voxel and the voxel spacing in mm. \n",
    "\n",
    "To find the voxel coordinates of a nodule, given its real world position, we use the ```GetOrigin()``` and ```GetSpacing()``` method of the itk image object:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "itk_img = sitk.ReadImage(img_file) \n",
    "img_array = sitk.GetArrayFromImage(itk_img) # indexes are z,y,x (notice the ordering)\n",
    "center = np.array([node_x,node_y,node_z])   # nodule center\n",
    "origin = np.array(itk_img.GetOrigin())      # x,y,z  Origin in world coordinates (mm)\n",
    "spacing = np.array(itk_img.GetSpacing())    # spacing of voxels in world coor. (mm)\n",
    "v_center =np.rint((center-origin)/spacing)  # nodule center in voxel space (still x,y,z ordering)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The center of the nodule is located in the ```v_center[2]``` slice of the img_array. We pass the node information to the ```make_mask()``` function and copy the generated masks and the image for the ```v_center[2]```slice and the slice above and below it. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "i = 0\n",
    "for i_z in range(int(v_center[2])-1,int(v_center[2])+2):\n",
    "    mask = make_mask(center,diam,i_z*spacing[2]+origin[2],width,height,spacing,origin)\n",
    "    masks[i] = mask\n",
    "    imgs[i] = matrix2int16(img_array[i_z])\n",
    "    i+=1\n",
    "np.save(output_path+\"images_%d.npy\" % (fcount) ,imgs)\n",
    "np.save(output_path+\"masks_%d.npy\" % (fcount) ,masks)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the ```make_mask()``` function it is worth noting that the mask coordinates have to match the ordering of the array coordinates. The ```x``` and ```y``` ordering is flipped. See the next to last line in the below code:\n",
    "\n",
    "```\n",
    "    def make_mask(center,diam,z,width,height,spacing,origin):\n",
    "    \n",
    "    ...\n",
    "    \n",
    "    for v_x in v_xrange:\n",
    "        for v_y in v_yrange:\n",
    "            p_x = spacing[0]*v_x + origin[0]\n",
    "            p_y = spacing[1]*v_y + origin[1]\n",
    "            if np.linalg.norm(center-np.array([p_x,p_y,z]))<=diam:\n",
    "                mask[int((p_y-origin[1])/spacing[1]),int((p_x-origin[0])/spacing[0])] = 1.0\n",
    "    return(mask)\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Should we collect more slices from each scan?\n",
    "\n",
    "Since the nodule locations are defined in terms of spheres, and the nodules are irregularly shaped, slices near the edges of the spheres may contain no nodule tissue. Using such slices would contaminate the training set with false positives. For this segmentation project, there is probably an optimal number of slices through a nodule that one should incorporate. For simplicity, we stick to 3 and only pull the slices centered on the largest nodule. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Check to make sure the nodule masks look as expected:\n",
    "\n",
    "```\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "imgs = np.load(output_path+'images_0.npy')\n",
    "masks = np.load(output_path+'masks_0.npy')\n",
    "for i in range(len(imgs)):\n",
    "    print \"image %d\" % i\n",
    "    fig,ax = plt.subplots(2,2,figsize=[8,8])\n",
    "    ax[0,0].imshow(imgs[i],cmap='gray')\n",
    "    ax[0,1].imshow(masks[i],cmap='gray')\n",
    "    ax[1,0].imshow(imgs[i]*masks[i],cmap='gray')\n",
    "    plt.show()\n",
    "    raw_input(\"hit enter to cont : \")\n",
    "```\n",
    "The image on the top left is the scan slice. The image on the top right is the node mask. The image on the bottom left is the masked slice, highlighting the node. \n",
    "\n",
    "![Example LUNA Mask](https://kaggle2.blob.core.windows.net/competitions/kaggle/6004/media/example_luna_mask.png)\n",
    "\n",
    "Close up on the nodule : \n",
    "\n",
    "![Example LUNA Mask detail](https://kaggle2.blob.core.windows.net/competitions/kaggle/6004/media/example_luna_mask_detail.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Isolation of the Lung Region of Interest to Narrow Our Nodule Search\n",
    "\n",
    "The node masks seem to be constructed properly. The next step is to isolate the lungs in the images. We'll need\n",
    "to import some skimage image processing modules for this step. The general strategy is to threshold the image to isolate the \n",
    "regions within the image, and then to identify which of those regions are the lungs. The lungs have a high constrast with the surrounding tissue, so the thresholding is fairly straightforward. We use some ad-hoc criteria for eliminating the non-lung regions from the image which do not apply equally to all data sets. \n",
    "\n",
    "In addition to our previous imports, we'll make use of ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "from skimage import morphology\n",
    "from skimage import measure\n",
    "from sklearn.cluster import KMeans\n",
    "from skimage.transform import resize\n",
    "```\n",
    "\n",
    "These steps are found in LUNA_segment_lung_ROI.py\n",
    "\n",
    "The arrays were loaded as ```dtype = np.float64``` because ```KMeans``` in sklearn has a bug related to the precision of the input to ```KMeans```. \n",
    "\n",
    "We'll walk through the steps of isolating the lung ROI with ```img```, which is a 512 X 512 slice from the set we extracted from the LUNA 2016 data, which looks like:\n",
    "\n",
    "![ROI Step 1](https://kaggle2.blob.core.windows.net/competitions/kaggle/6004/media/ROI_step1.png)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Thresholding\n",
    "\n",
    "Our first step is to standardize the pixel values and take a look at the intensity distribution\n",
    "\n",
    "```\n",
    "img = imgs_to_process[i]\n",
    "#Standardize the pixel values\n",
    "mean = np.mean(img)\n",
    "std = np.std(img)\n",
    "img = img-mean\n",
    "img = img/std\n",
    "plt.hist(img.flatten(),bins=200)\n",
    "```\n",
    "![ROI Step 1 hist](https://kaggle2.blob.core.windows.net/competitions/kaggle/6004/media/ROI_step1_hist.png)\n",
    "\n",
    "The underflow peak near -1.5 is the black out-of-scanner part of the image. The peaks around 0.0 are the background and lung interior and the wide clumps from 1.0 to 2.0 are the non-lung-tissue and bone. The structure of this histogram varies throughout the data set. Two images are shown below that are typical of the data set. The one on the left has the same black background around a grey circular region of scanner data as is present in ```img```. That black background is not present in the image on the right, making for a very different pixel value histogram.   \n",
    "\n",
    "![ROI Hist Diff](https://kaggle2.blob.core.windows.net/competitions/kaggle/6004/media/ROI_hist_diff.png)\n",
    "\n",
    "We have to make sure that we set our threshold between the lung pixel values and the denser tissue pixel values. To do this, we reset the pixels with the minimum value to the average pixel value near the center of the picture and perform kmeans clustering with k=2. This seems to work well for both scenarios. \n",
    "\n",
    "```\n",
    "middle = img[100:400,100:400] \n",
    "mean = np.mean(middle)  \n",
    "max = np.max(img)\n",
    "min = np.min(img)\n",
    "#move the underflow bins\n",
    "img[img==max]=mean\n",
    "img[img==min]=mean\n",
    "\n",
    "kmeans = KMeans(n_clusters=2).fit(np.reshape(middle,[np.prod(middle.shape),1]))\n",
    "centers = sorted(kmeans.cluster_centers_.flatten())\n",
    "threshold = np.mean(centers)\n",
    "thresh_img = np.where(img<threshold,1.0,0.0)  # threshold the image\n",
    "```\n",
    "\n",
    "Which produces a satisfactory separation of regions for both types of images and eliminates the black halo in the one on the left\n",
    "\n",
    "![ROI Step 2](https://kaggle2.blob.core.windows.net/competitions/kaggle/6004/media/ROI_step2.png)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Erosion and Dilation\n",
    "\n",
    "We then use an erosion and dilation to fill in the incursions into the lungs region by radio-opaque tissue, followed by a selection of the regions based on the bounding box sizes of each region. The initial set of regions looks like \n",
    "```\n",
    "eroded = morphology.erosion(thresh_img,np.ones([4,4]))\n",
    "dilation = morphology.dilation(eroded,np.ones([10,10]))\n",
    "labels = measure.label(dilation)\n",
    "label_vals = np.unique(labels)\n",
    "plt.imshow(labels)\n",
    "```\n",
    "![ROI Step 3](https://kaggle2.blob.core.windows.net/competitions/kaggle/6004/media/ROI_step3.png)\n",
    "\n",
    "### Cutting non-ROI Regions\n",
    "The cuts applied to each region bounding box were determined empirically and seem to work well for the LUNA data, but may not be generally applicable\n",
    "```\n",
    "labels = measure.label(dilation)\n",
    "label_vals = np.unique(labels)\n",
    "regions = measure.regionprops(labels)\n",
    "good_labels = []\n",
    "for prop in regions:\n",
    "    B = prop.bbox\n",
    "    if B[2]-B[0]<475 and B[3]-B[1]<475 and B[0]>40 and B[2]<472:\n",
    "        good_labels.append(prop.label)\n",
    "mask = np.ndarray([512,512],dtype=np.int8)\n",
    "mask[:] = 0\n",
    "#\n",
    "#  The mask here is the mask for the lungs--not the nodes\n",
    "#  After just the lungs are left, we do another large dilation\n",
    "#  in order to fill in and out the lung mask \n",
    "#\n",
    "for N in good_labels:\n",
    "    mask = mask + np.where(labels==N,1,0)\n",
    "mask = morphology.dilation(mask,np.ones([10,10])) # one last dilation\n",
    "plt.imshow(mask,cmap='gray')\n",
    "```\n",
    "![ROI Step 4](https://kaggle2.blob.core.windows.net/competitions/kaggle/6004/media/ROI_step4.png)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Applying the ROI Masks\n",
    "\n",
    "The next step in LUNA_segment_lung_ROI.py is applying the mask of the lung ROI to each of the images, cropping down to the bounding square of the lungs ROI, and then resizing the resulting image to 512 X 512. \n",
    "```\n",
    "masks = np.load(working_path+\"lungmask_0.py\")\n",
    "imgs = np.load(working_path+\"images_0.py\")\n",
    "imgs = masks*imgs\n",
    "```\n",
    "...crop to bounding square and resize to 512 X 512...\n",
    "\n",
    "Then we perform some final pixel normalization. This is because the mask sends the non ROI area in the picture to 0, and that operation is not sensitive to the pixel value distribution. To fix this, we get the mean and standard deviation of the masked region and send the background (now zero) to the lower end of the pixel distribution (```-1.2*stdev```, which was chosen empirically). \n",
    "```\n",
    "#\n",
    "# renormalizing the masked image (in the mask region)\n",
    "#\n",
    "new_mean = np.mean(img[mask>0])  \n",
    "new_std = np.std(img[mask>0])\n",
    "#\n",
    "#  Pushing the background color up to the lower end\n",
    "#  of the pixel range for the lungs\n",
    "#\n",
    "old_min = np.min(img)       # background color\n",
    "img[img==old_min] = new_mean-1.2*new_std   # resetting backgound color\n",
    "img = img-new_mean\n",
    "img = img/new_std\n",
    "```\n",
    "\n",
    "The final product is a set of lungs that is ready to be compiled into our training example set. \n",
    "\n",
    "![ROI Step 5](https://kaggle2.blob.core.windows.net/competitions/kaggle/6004/media/ROI_step5.png)\n",
    "\n",
    "These images and the correspondingly trimmed and rescaled masks are randomized and sent to a single file that contains a numpy array of dimension ```[<num_images>,1,512,512]```. The 1 is important as the U-net is enabled for multiple channels. \n",
    "\n",
    "```\n",
    "#\n",
    "#  Writing out images and masks as 1 channel arrays for input into network\n",
    "#\n",
    "final_images = np.ndarray([num_images,1,512,512],dtype=np.float32)\n",
    "final_masks = np.ndarray([num_images,1,512,512],dtype=np.float32)\n",
    "for i in range(num_images):\n",
    "    final_images[i,0] = out_images[i]\n",
    "    final_masks[i,0] = out_nodemasks[i]\n",
    "\n",
    "rand_i = np.random.choice(range(num_images),size=num_images,replace=False)\n",
    "test_i = int(0.2*num_images)\n",
    "np.save(working_path+\"trainImages.npy\",final_images[rand_i[test_i:]])\n",
    "np.save(working_path+\"trainMasks.npy\",final_masks[rand_i[test_i:]])\n",
    "np.save(working_path+\"testImages.npy\",final_images[rand_i[:test_i]])\n",
    "np.save(working_path+\"testMasks.npy\",final_masks[rand_i[:test_i]])\n",
    "\n",
    "```\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can check the ROI isolationg by examining the lung masks and the original files:\n",
    "\n",
    "```\n",
    "imgs = np.load(working path+'images_0.npy')\n",
    "lungmask = np.load(working_path+'lungmask_0.npy')\n",
    "\n",
    "for i in range(len(imgs)):\n",
    "    print \"image %d\" % i\n",
    "    fig,ax = plt.subplots(2,2,figsize=[8,8])\n",
    "    ax[0,0].imshow(imgs[i],cmap='gray')\n",
    "    ax[0,1].imshow(lungmask[i],cmap='gray')\n",
    "    ax[1,0].imshow(imgs[i]*lungmask[i],cmap='gray')\n",
    "    plt.show()\n",
    "    raw_input(\"hit enter to cont : \")\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dice Ceofficient Cost function for Segmentation\n",
    "The network we'll be using is the U-net linked to the beginning of the tutorial which \n",
    "uses the keras framework. The loss function is the dice coefficient \n",
    "https://en.wikipedia.org/wiki/S%C3%B8rensen%E2%80%93Dice_coefficient\n",
    "comparing the predicted and actual node mask. \n",
    "\n",
    "The following code snippets are all taken from LUNA_train_unet.py\n",
    "\n",
    "The loss function is as follows: \n",
    "\n",
    "```\n",
    "smooth = 1.\n",
    "# Tensorflow version for the model\n",
    "def dice_coef(y_true, y_pred):\n",
    "    y_true_f = K.flatten(y_true)\n",
    "    y_pred_f = K.flatten(y_pred)\n",
    "    intersection = K.sum(y_true_f * y_pred_f)\n",
    "    return (2. * intersection + smooth) / (K.sum(y_true_f) + K.sum(y_pred_f) + smooth)\n",
    "def dice_coef_loss(y_true, y_pred):\n",
    "    return -dice_coef(y_true, y_pred)\n",
    "```\n",
    "Which is similar to the metric used to evaluate the Ultrasound Nerve Segmentation challenge that this network was originally written for (once again, see the link at the beginning of this tutorial).\n",
    "\n",
    "### Loading the Segmenter\n",
    "The function call\n",
    "```\n",
    "model = get_unet() \n",
    "model_checkpoint = ModelCheckpoint('unet.hdf5', monitor='loss', save_best_only=True)\n",
    "```\n",
    "will compile and return the model and tell keras to save the model weights during checkpoints. If you want to load the best weights from a previous training session or use the weights included in this tutorial's repo, load the weight file with the line\n",
    "```\n",
    "model.load_weights('unet.hdf5')\n",
    "```\n",
    "\n",
    "### Training the Segmenter\n",
    "Calling LUNA_train_unet.py from the command line will automatically attempt to load a unet.hdf5 file from the current directory, train it according to the parameters set by the line\n",
    "```\n",
    "model.fit(imgs_train, imgs_mask_train, batch_size=2, nb_epoch=20, \n",
    "           verbose=1, shuffle=True,callbacks=[model_checkpoint])\n",
    "```\n",
    "in the script, and make predictions on the test set.\n",
    "```\n",
    "num_test = len(imgs_test)\n",
    "imgs_mask_test = np.ndarray([num_test,1,512,512],dtype=np.float32)\n",
    "for i in range(num_test):\n",
    "    imgs_mask_test[i] = model.predict([imgs_test[i:i+1]], verbose=0)[0]\n",
    "np.save('masksTestPredicted.npy', imgs_mask_test)\n",
    "mean = 0.0\n",
    "for i in range(num_test):\n",
    "    mean+=dice_coef_np(imgs_mask_test_true[i,0], imgs_mask_test[i,0])\n",
    "mean/=num_test\n",
    "print(\"Mean Dice Coeff : \",mean)\n",
    "```\n",
    "The ```model.predict()``` function can take more than one case at a time, but that can quickly overload a GPU, so we are looping over individual cases. \n",
    "\n",
    "The final results for this tutorial were produced using a multi-GPU machine using TitanX's. For a home GPU computation benchmark, a personal set up with a GTX970 we were able to run 20 epochs with a training set size of 320 and batch size of 2 in about an hour. We started obtaing reasonable nodule mask predictions after about 3 hours of training once the reported loss value approached 0.3. \n",
    "\n",
    "An example segmentation is given here for the three slices taken from a patient scan. The perfect circle is the \"true\" node mask from the LUNA annotationc.csv file, and the red is the predicted node region from the segmenter. The original image is given in the top right. \n",
    "\n",
    "![Example Segmentation](https://kaggle2.blob.core.windows.net/competitions/kaggle/6004/media/example_segmentation.png)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Training a Classifier for Identifying Cancer\n",
    "\n",
    "Now we are ready to begin training a classifier using our image segmentation from the previous sections to generate features. \n",
    "\n",
    "The Data Science Bowl training data set must be fed through the segmenter, which can be done largely by reusing the code used to treat the LUNA data. There are two points where this deviates.\n",
    "\n",
    "First of all, the DSB data is in dicom format, which can be read using the pydicom module\n",
    "```\n",
    "import dicom\n",
    "dc = dicom.read_file(filename)\n",
    "img = dc.pixel_array\n",
    "```\n",
    "Secondly, in order to locate nodes in the scans, every layer of the scan must be run through the segmenter, and thus every layer must also be subject to the image processing to mask off the ROI. This can be a very time consuming process. \n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple classifier based on nodule features\n",
    "We start by characterizing some of the features of the nodule maps and putting them into a feature vector that we can use for classification purposes. Our list of features is by no means exhaustive and is meant to illustrate the process of developing some meterics that could be used to characterize the segmented regions where nodules are likely.\n",
    "\n",
    "We encourage you to play with adding some new features, explore convolutional models for extracting features from the region of interest directly. We have included some features about the average size, morphology, and position within the image for use in model building. \n",
    "\n",
    "\n",
    "```\n",
    "def getRegionMetricRow(fname = \"nodules.npy\"):\n",
    "    seg = np.load(fname)\n",
    "    nslices = seg.shape[0]\n",
    "    \n",
    "    #metrics\n",
    "    totalArea = 0.\n",
    "    avgArea = 0.\n",
    "    maxArea = 0.\n",
    "    avgEcc = 0.\n",
    "    avgEquivlentDiameter = 0.\n",
    "    stdEquivlentDiameter = 0.\n",
    "    weightedX = 0.\n",
    "    weightedY = 0.\n",
    "    numNodes = 0.\n",
    "    numNodesperSlice = 0.\n",
    "    # do not allow any nodes to be larger than 10% of the pixels to eliminate background regions\n",
    "    maxAllowedArea = 0.10 * 512 * 512 \n",
    "    \n",
    "    areas = []\n",
    "    eqDiameters = []\n",
    "    for slicen in range(nslices):\n",
    "        regions = getRegionFromMap(seg[slicen,0,:,:])\n",
    "        for region in regions:\n",
    "            if region.area > maxAllowedArea:\n",
    "                continue\n",
    "            totalArea += region.area\n",
    "            areas.append(region.area)\n",
    "            avgEcc += region.eccentricity\n",
    "            avgEquivlentDiameter += region.equivalent_diameter\n",
    "            eqDiameters.append(region.equivalent_diameter)\n",
    "            weightedX += region.centroid[0]*region.area\n",
    "            weightedY += region.centroid[1]*region.area\n",
    "            numNodes += 1\n",
    "            \n",
    "    weightedX = weightedX / totalArea \n",
    "    weightedY = weightedY / totalArea\n",
    "    avgArea = totalArea / numNodes\n",
    "    avgEcc = avgEcc / numNodes\n",
    "    avgEquivlentDiameter = avgEquivlentDiameter / numNodes\n",
    "    stdEquivlentDiameter = np.std(eqDiameters)\n",
    "    \n",
    "    maxArea = max(areas)\n",
    "    \n",
    "    \n",
    "    numNodesperSlice = numNodes*1. / nslices\n",
    "    \n",
    "    \n",
    "    return np.array([avgArea,maxArea,avgEcc,avgEquivlentDiameter,\\\n",
    "                     stdEquivlentDiameter, weightedX, weightedY, numNodes, numNodesperSlice])\n",
    "```\n",
    "\n",
    "```\n",
    "def getRegionFromMap(slice_npy):\n",
    "    thr = np.where(slice_npy > np.mean(slice_npy),0.,1.0)\n",
    "    label_image = label(thr)\n",
    "    labels = label_image.astype(int)\n",
    "    regions = regionprops(labels)\n",
    "    return regions\n",
    "```\n",
    "\n",
    "```\n",
    "import pickle\n",
    "def createFeatureDataset(nodfiles=None):\n",
    "    if nodfiles == None:\n",
    "        noddir = \"/home/jmulholland/NLST_nodules/\"\n",
    "        nodfiles = glob(noddir +\"*npy\")\n",
    "    # dict with mapping between truth and \n",
    "    truthdata = pickle.load(open(\"/home/sander/truthdict.pkl\",'r'))\n",
    "    numfeatures = 9\n",
    "    feature_array = np.zeros((len(nodfiles),numfeatures))\n",
    "    truth_metric = np.zeros((len(nodfiles)))\n",
    "    \n",
    "    for i,nodfile in enumerate(nodfiles):\n",
    "        patID = nodfile.split(\"_\")[2]\n",
    "        truth_metric[i] = truthdata[int(patID)]\n",
    "        feature_array[i] = getRegionMetricRow(nodfile)\n",
    "    \n",
    "    np.save(\"dataY.npy\", truth_metric)\n",
    "    np.save(\"dataX.npy\", feature_array)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we've created feature vectors, we load them into some simple classification models and see how we're performing. We've choosen random forest and XGBoost to create some models based on our feature engineering. \n",
    "\n",
    "```\n",
    "from sklearn import cross_validation\n",
    "from sklearn.cross_validation import StratifiedKFold as KFold\n",
    "from sklearn.metrics import classification_report\n",
    "from sklearn.ensemble import RandomForestClassifier as RF\n",
    "import xgboost as xgb\n",
    "\n",
    "X = np.load(\"dataX.npy\")\n",
    "Y = np.load(\"dataY.npy\")\n",
    "\n",
    "kf = KFold(Y, n_folds=3)\n",
    "y_pred = Y * 0\n",
    "for train, test in kf:\n",
    "    X_train, X_test, y_train, y_test = X[train,:], X[test,:], Y[train], Y[test]\n",
    "    clf = RF(n_estimators=100, n_jobs=3)\n",
    "    clf.fit(X_train, y_train)\n",
    "    y_pred[test] = clf.predict(X_test)\n",
    "print classification_report(Y, y_pred, target_names=[\"No Cancer\", \"Cancer\"])\n",
    "print(\"logloss\",logloss(Y, y_pred))\n",
    "\n",
    "# All Cancer\n",
    "print \"Predicting all positive\"\n",
    "y_pred = np.ones(Y.shape)\n",
    "print classification_report(Y, y_pred, target_names=[\"No Cancer\", \"Cancer\"])\n",
    "print(\"logloss\",logloss(Y, y_pred))\n",
    "\n",
    "# No Cancer\n",
    "print \"Predicting all negative\"\n",
    "y_pred = Y*0\n",
    "print classification_report(Y, y_pred, target_names=[\"No Cancer\", \"Cancer\"])\n",
    "print(\"logloss\",logloss(Y, y_pred))\n",
    "\n",
    "# try XGBoost\n",
    "print (\"XGBoost\")\n",
    "kf = KFold(Y, n_folds=3)\n",
    "y_pred = Y * 0\n",
    "for train, test in kf:\n",
    "    X_train, X_test, y_train, y_test = X[train,:], X[test,:], Y[train], Y[test]\n",
    "    clf = xgb.XGBClassifier(objective=\"binary:logistic\")\n",
    "    clf.fit(X_train, y_train)\n",
    "    y_pred[test] = clf.predict(X_test)\n",
    "print classification_report(Y, y_pred, target_names=[\"No Cancer\", \"Cancer\"])\n",
    "print(\"logloss\",logloss(Y, y_pred))\n",
    "\n",
    "```\n",
    "The output shows the results for these modeling efforts. \n",
    "\n",
    "```Random Forest\n",
    "            precision    recall  f1-score   support\n",
    "\n",
    "  No Cancer       0.77      0.95      0.85       219\n",
    "     Cancer       0.08      0.02      0.03        64\n",
    "\n",
    "avg / total       0.61      0.73      0.66       283\n",
    "\n",
    "('logloss', 9.1534198755740697)\n",
    "Predicting all positive\n",
    "             precision    recall  f1-score   support\n",
    "\n",
    "  No Cancer       0.00      0.00      0.00       219\n",
    "     Cancer       0.23      1.00      0.37        64\n",
    "\n",
    "avg / total       0.05      0.23      0.08       283\n",
    "\n",
    "('logloss', 26.728505803260379)\n",
    "Predicting all negative\n",
    "             precision    recall  f1-score   support\n",
    "\n",
    "  No Cancer       0.77      1.00      0.87       219\n",
    "     Cancer       0.00      0.00      0.00        64\n",
    "\n",
    "avg / total       0.60      0.77      0.68       283\n",
    "\n",
    "('logloss', 7.8108893613932242)\n",
    "XGBoost\n",
    "             precision    recall  f1-score   support\n",
    "\n",
    "  No Cancer       0.77      0.92      0.84       219\n",
    "     Cancer       0.18      0.06      0.09        64\n",
    "\n",
    "avg / total       0.64      0.72      0.67       283\n",
    "\n",
    "('logloss', 9.5195722669850706)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Once we've created our feature dataset, we can create a simple classifer based on the featureset to determine the cancer diagnosis. We've chose to illustrate two popular classifiers Random Forest and XGBoost, because they are both robust to overfitting and quick to train. The common python implementation of XGBoost also allows for easy modification of the objective function and reweighting the class importance, so that we can reweight the optimization problem to give more weight to cancer cases.\n",
    "\n",
    "```\n",
    "from sklearn import cross_validation\n",
    "from sklearn.cross_validation import StratifiedKFold as KFold\n",
    "from sklearn.metrics import classification_report\n",
    "from sklearn.ensemble import RandomForestClassifier as RF\n",
    "import xgboost as xgb\n",
    "\n",
    "X = np.load(\"dataX.npy\")\n",
    "Y = np.load(\"dataY.npy\")\n",
    "\n",
    "kf = KFold(Y, n_folds=3)\n",
    "y_pred = Y * 0\n",
    "for train, test in kf:\n",
    "    X_train, X_test, y_train, y_test = X[train,:], X[test,:], Y[train], Y[test]\n",
    "    clf = RF(n_estimators=100, n_jobs=3)\n",
    "    clf.fit(X_train, y_train)\n",
    "    y_pred[test] = clf.predict(X_test)\n",
    "print(\"Random Forest\")\n",
    "print classification_report(Y, y_pred, target_names=[\"No Cancer\", \"Cancer\"])\n",
    "print(\"logloss\",logloss(Y, y_pred))\n",
    "\n",
    "# try XGBoost\n",
    "print (\"XGBoost\")\n",
    "kf = KFold(Y, n_folds=3)\n",
    "y_pred = Y * 0\n",
    "for train, test in kf:\n",
    "    X_train, X_test, y_train, y_test = X[train,:], X[test,:], Y[train], Y[test]\n",
    "    clf = xgb.XGBClassifier(objective=\"binary:logistic\", scale_pos_weight=3)\n",
    "    clf.fit(X_train, y_train)\n",
    "    y_pred[test] = clf.predict(X_test)\n",
    "print classification_report(Y, y_pred, target_names=[\"No Cancer\", \"Cancer\"])\n",
    "print(\"logloss\",logloss(Y, y_pred))\n",
    "\n",
    "# All Cancer\n",
    "print \"Predicting all positive\"\n",
    "y_pred = np.ones(Y.shape)\n",
    "print classification_report(Y, y_pred, target_names=[\"No Cancer\", \"Cancer\"])\n",
    "print(\"logloss\",logloss(Y, y_pred))\n",
    "\n",
    "# No Cancer\n",
    "print \"Predicting all negative\"\n",
    "y_pred = Y*0\n",
    "print classification_report(Y, y_pred, target_names=[\"No Cancer\", \"Cancer\"])\n",
    "print(\"logloss\",logloss(Y, y_pred))\n",
    "\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compare the results for Random forest, XGBoost, and two models consisting of only predicting cancer and only predicting no cancer. \n",
    "\n",
    "Output:\n",
    "```\n",
    "Random Forest\n",
    "             precision    recall  f1-score   support\n",
    "\n",
    "  No Cancer       0.81      0.97      0.88       463\n",
    "     Cancer       0.13      0.02      0.03       107\n",
    "\n",
    "avg / total       0.68      0.79      0.72       570\n",
    "\n",
    "('logloss', 7.150150893624649)\n",
    "XGBoost\n",
    "             precision    recall  f1-score   support\n",
    "\n",
    "  No Cancer       0.83      0.86      0.84       463\n",
    "     Cancer       0.27      0.21      0.24       107\n",
    "\n",
    "avg / total       0.72      0.74      0.73       570\n",
    "\n",
    "('logloss', 8.9074570257718957)\n",
    "\n",
    "Predicting all positive\n",
    "             precision    recall  f1-score   support\n",
    "\n",
    "  No Cancer       0.00      0.00      0.00       463\n",
    "     Cancer       0.19      1.00      0.32       107\n",
    "\n",
    "avg / total       0.04      0.19      0.06       570\n",
    "\n",
    "('logloss', 28.055831025357818)\n",
    "Predicting all negative\n",
    "             precision    recall  f1-score   support\n",
    "\n",
    "  No Cancer       0.81      1.00      0.90       463\n",
    "     Cancer       0.00      0.00      0.00       107\n",
    "\n",
    "avg / total       0.66      0.81      0.73       570\n",
    "\n",
    "('logloss', 6.4835948671148085)\n",
    "\n",
    "\n",
    "```\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Where to go next?\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "We've given you a framework for approaching this problem that combines a deep learning based segmentation approach with an older computer vision approach of hand designed features. From here you have lots of ways to improve the u-net model with more data, or additional preprocessing. The classification piece could be replaced by a another convolutional nueral net or you could take more advantage of the 3d nature of the nodules that we've mostly considered as independant 2d slices here."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}