{ "cells": [ { "cell_type": "markdown", "metadata": { "_cell_guid": "efd5feb5-8262-1b45-af15-0dd16116052b" }, "source": [ "In this kernel, I will be talking about the methods that will help in better understanding of the problem statement and visualisation of the data. I will also provide links of useful resources and information.\n", "\n", "The script is written in Python. I recommend people to install anaconda on their desktop because of its advantages mentioned [here][1]. The libraries that will be used in this tutorial for reading, processing and visualisation of data are matplotlib, numpy, skimage and [pydicom][2]. \n", "\n", "The images are of size (z, 512, 512) where z is the number of slices in the CT Scan and varies depending on the resolution of the scanner. Such large images cant be fed directly into a Convolution Network architectures because of the limit on the computation power. Thus we will have to find the regions that are more probable of having cancer. We will be reducing our search space by first segmenting the lungs and then removing the low intensity regions. \n", "\n", "In this tutorial, we will first start with reading the dataset and visualising it. After that, we will be segmenting the lung structures and then find the region of interest(possible cancer regions) in the CT Scans using Image processing methods. Then I will talk about how to preprocess LUNA16 dataset for training architectures like UNet for segmentation and candidate classification.\n", "\n", "The segmentation of lung structures is very challenging problem because homogeneity is not present in the lung region, similar densities in the pulmonary structures, different scanners and scanning protocols. The segmented lungs can be further used to find the lung nodule candidates and regions of interest which may help in better classification of the CT Scans. Finding the lung nodule regions is a very hard problem because there are nodules that are attached to the blood vessels or are present at the boundary of the lung region. The lung nodule candidates can be further used for classification by cutting 3D voxels around them and passing it through a 3D CNNs which can be trained on LUNA16 dataset. The LUNA 16 dataset has the location of the nodules in each CT Scan thus will be useful for training the classifier. \n", "\n", " [1]: https://www.reddit.com/r/Python/comments/3t23vv/what_advantages_are_there_of_using_anaconda/\n", " [2]: http://pydicom.readthedocs.io/en/stable/getting_started.html" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "ceec3ade-181c-4059-e0f2-958a1f7629af" }, "source": [ "**Reading a CT Scan**\n", "-----------------\n", "The input folder has three things, one is the sample_images folders which has the sample CT Scans. The `stage1_labels.csv` contains the cancer ground truth for the stage 1 training set images and `stage1_sample_submission.csv` shows the submission format for stage 1. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "53553842-54eb-d04d-1f5a-fd0c08472863" }, "outputs": [], "source": [ "# This Python 3 environment comes with many helpful analytics libraries installed\n", "# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n", "# For example, here's several helpful packages to load in \n", "\n", "import numpy as np # linear algebra\n", "import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n", "import skimage, os\n", "from skimage.morphology import ball, disk, dilation, binary_erosion, remove_small_objects, erosion, closing, reconstruction, binary_closing\n", "from skimage.measure import label,regionprops, perimeter\n", "from skimage.morphology import binary_dilation, binary_opening\n", "from skimage.filters import roberts, sobel\n", "from skimage import measure, feature\n", "from skimage.segmentation import clear_border\n", "from skimage import data\n", "from scipy import ndimage as ndi\n", "import matplotlib.pyplot as plt\n", "from mpl_toolkits.mplot3d.art3d import Poly3DCollection\n", "import dicom\n", "import scipy.misc\n", "import numpy as np\n", "\n", "# Input data files are available in the \"../input/\" directory.\n", "# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n", "\n", "from subprocess import check_output\n", "print(check_output([\"ls\", \"../input\"]).decode(\"utf8\"))" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "13bdd5df-5f74-f085-e590-57c135ccd704" }, "source": [ "Each 3D CT Scan consists of many slices, whose number depends on the resolution of the scanner and each slice has a Instance Number associated with it which tells the index of the slice from the top. All the dicom files for a CT Scan are inside one folder having the CT Scan's name. Now we will read all the dicom slices for a scan and then stack them with respect to their Instance Number to get the 3D Lung CT Scanned Image." ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "0c33fe52-ca9f-477a-a85f-04d146963fc2" }, "source": [ "The sample_images folder has around 20 folders each corresponding to one CT Scan. Inside the folder there are many dicom files. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "9223bd08-30c1-6841-e4af-9be561bf756a" }, "outputs": [], "source": [ "print(check_output([\"ls\", \"../input/sample_images/\"]).decode(\"utf8\"))" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "bfe2dd1f-6061-b2dd-94a5-559a70d09f5b" }, "source": [ "Each CT Scan consists of multiple 2D slices which are provided in a DICOM format. At first, I will read the random dicom file of a CT Scan. After reading the image file, we will update the intensity values of -2000 with 0 because they are the pixels that fall outside of the scanner bounds." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "a4a6a3e0-5766-9421-5a9f-d9237773f860" }, "outputs": [], "source": [ "# Any results you write to the current directory are saved as output.\n", "lung = dicom.read_file('../input/sample_images/00cba091fa4ad62cc3200a657aeb957e/38c4ff5d36b5a6b6dc025435d62a143d.dcm')\n", "\n", "slice = lung.pixel_array\n", "slice[slice == -2000] = 0\n", "plt.imshow(slice, cmap=plt.cm.gray)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "740a14d6-edc1-cfbd-c31c-6a378a5d71be" }, "outputs": [], "source": [ "def read_ct_scan(folder_name):\n", " # Read the slices from the dicom file\n", " slices = [dicom.read_file(folder_name + filename) for filename in os.listdir(folder_name)]\n", " \n", " # Sort the dicom slices in their respective order\n", " slices.sort(key=lambda x: int(x.InstanceNumber))\n", " \n", " # Get the pixel values for all the slices\n", " slices = np.stack([s.pixel_array for s in slices])\n", " slices[slices == -2000] = 0\n", " return slices" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "9b923da3-6985-27f1-0cef-f8c6dd0434ad" }, "outputs": [], "source": [ "ct_scan = read_ct_scan('../input/sample_images/00cba091fa4ad62cc3200a657aeb957e/') " ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "8db1af50-854f-da03-3c03-7782b7fd39d9" }, "source": [ "To visualise the slices, we will have to plot them. `matplotlib` is used for plotting the slices. The `plot_ct_scan` function takes a 3D CT Scanned Image array as input and plots equally spaced slices. The CT Scans are grayscale images i.e. the value of each pixel is a single sample, which means it carries only intensity information." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "43b4d7bc-66d8-b12a-47e4-d391e998f796" }, "outputs": [], "source": [ "def plot_ct_scan(scan):\n", " f, plots = plt.subplots(int(scan.shape[0] / 20) + 1, 4, figsize=(25, 25))\n", " for i in range(0, scan.shape[0], 5):\n", " plots[int(i / 20), int((i % 20) / 5)].axis('off')\n", " plots[int(i / 20), int((i % 20) / 5)].imshow(scan[i], cmap=plt.cm.bone) " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "0f3389c1-65b9-8af3-edc8-6aa38b1e6498" }, "outputs": [], "source": [ "plot_ct_scan(ct_scan)" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "4118db22-f1d7-304a-a58f-0f8aaca8311a" }, "source": [ "## Segmentation of Lungs ##\n", "After reading the CT Scan, the first step in preprocessing is the segmentation of lung structures because it is obvious that the regions of interests lies inside the lungs. It is visible that the lungs are the darker regions in the CT Scans. The bright region inside the lungs are the blood vessels or air. A threshold of 604(-400 HU) is used at all places because it was found in experiments that it works just fine. We segment lung structures from each slice of the CT Scan image and try not to loose the possible region of interests attached to the lung wall. There are some nodules which may be attached to the lung wall.\n", "\n", "I will first explain a common method using simple Image Processing and Morphological operations to segment the lungs and then will give references and summaries to good links of papers. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "3fc2eca2-e580-e66e-561b-d1812e829897" }, "outputs": [], "source": [ "def get_segmented_lungs(im, plot=False):\n", " \n", " '''\n", " This funtion segments the lungs from the given 2D slice.\n", " '''\n", " if plot == True:\n", " f, plots = plt.subplots(8, 1, figsize=(5, 40))\n", " '''\n", " Step 1: Convert into a binary image. \n", " '''\n", " binary = im < 604\n", " if plot == True:\n", " plots[0].axis('off')\n", " plots[0].imshow(binary, cmap=plt.cm.bone) \n", " '''\n", " Step 2: Remove the blobs connected to the border of the image.\n", " '''\n", " cleared = clear_border(binary)\n", " if plot == True:\n", " plots[1].axis('off')\n", " plots[1].imshow(cleared, cmap=plt.cm.bone) \n", " '''\n", " Step 3: Label the image.\n", " '''\n", " label_image = label(cleared)\n", " if plot == True:\n", " plots[2].axis('off')\n", " plots[2].imshow(label_image, cmap=plt.cm.bone) \n", " '''\n", " Step 4: Keep the labels with 2 largest areas.\n", " '''\n", " areas = [r.area for r in regionprops(label_image)]\n", " areas.sort()\n", " if len(areas) > 2:\n", " for region in regionprops(label_image):\n", " if region.area < areas[-2]:\n", " for coordinates in region.coords: \n", " label_image[coordinates[0], coordinates[1]] = 0\n", " binary = label_image > 0\n", " if plot == True:\n", " plots[3].axis('off')\n", " plots[3].imshow(binary, cmap=plt.cm.bone) \n", " '''\n", " Step 5: Erosion operation with a disk of radius 2. This operation is \n", " seperate the lung nodules attached to the blood vessels.\n", " '''\n", " selem = disk(2)\n", " binary = binary_erosion(binary, selem)\n", " if plot == True:\n", " plots[4].axis('off')\n", " plots[4].imshow(binary, cmap=plt.cm.bone) \n", " '''\n", " Step 6: Closure operation with a disk of radius 10. This operation is \n", " to keep nodules attached to the lung wall.\n", " '''\n", " selem = disk(10)\n", " binary = binary_closing(binary, selem)\n", " if plot == True:\n", " plots[5].axis('off')\n", " plots[5].imshow(binary, cmap=plt.cm.bone) \n", " '''\n", " Step 7: Fill in the small holes inside the binary mask of lungs.\n", " '''\n", " edges = roberts(binary)\n", " binary = ndi.binary_fill_holes(edges)\n", " if plot == True:\n", " plots[6].axis('off')\n", " plots[6].imshow(binary, cmap=plt.cm.bone) \n", " '''\n", " Step 8: Superimpose the binary mask on the input image.\n", " '''\n", " get_high_vals = binary == 0\n", " im[get_high_vals] = 0\n", " if plot == True:\n", " plots[7].axis('off')\n", " plots[7].imshow(im, cmap=plt.cm.bone) \n", " \n", " return im" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "8d027175-ccd6-aa1d-ca26-fa70857ea9aa" }, "source": [ "The `get_segmented_lungs` function segments a 2D slice of the CT Scan. I have outputted the slice after all steps for better visualisation and understanding of the code and applied operations." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "0b95efee-2a29-f0c6-3728-4f4425123573" }, "outputs": [], "source": [ "get_segmented_lungs(ct_scan[71], True)" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "ffb87321-c5fc-c04d-2559-b2b2b27d9066" }, "source": [ "Now, I will segment the whole CT Scan slice by slice and show some slices of the CT Scan." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "b3048b12-2ee0-d144-a742-014998bb0458" }, "outputs": [], "source": [ "def segment_lung_from_ct_scan(ct_scan):\n", " return np.asarray([get_segmented_lungs(slice) for slice in ct_scan])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "75f42d28-320f-eeb9-5347-9049e8dcf85c" }, "outputs": [], "source": [ "segmented_ct_scan = segment_lung_from_ct_scan(ct_scan)\n", "plot_ct_scan(segmented_ct_scan)" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "10c687a9-2eec-9cb3-94e9-48717e1d99f2" }, "source": [ "Nodule Candidate/Region of Interest Generation\n", "---------------------------\n", "After segmenting the lung structures from the CT Scanned images, our task is to find the candidate regions with nodules since the search space is very large. Also, whole image can't be classified directly using 3D CNNs due to limit on computation, we need to find possible regions of cancer and then classify them. It was found in experiments that all the region of interests have intensity > 604(-400 HU). So, we used this threshold to filter the darker regions. This reduces the number of candidates by a large number and preserves all the important regions with high recall. We then classify all the candidate points to reduce the False Positives." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "5d890715-8444-a56a-5c23-9dbe8a2353a6" }, "outputs": [], "source": [ "segmented_ct_scan[segmented_ct_scan < 604] = 0\n", "plot_ct_scan(segmented_ct_scan)" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "8d8e510b-f23d-aeab-3fb0-76f14037942a" }, "source": [ "After filtering, there are still lot of noise because of blood vessels. Thus we further remove the two largest connected component." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "8d7c6201-1143-ec15-79da-a4d30f74d2da" }, "outputs": [], "source": [ "selem = ball(2)\n", "binary = binary_closing(segmented_ct_scan, selem)\n", "\n", "label_scan = label(binary)\n", "\n", "areas = [r.area for r in regionprops(label_scan)]\n", "areas.sort()\n", "\n", "for r in regionprops(label_scan):\n", " max_x, max_y, max_z = 0, 0, 0\n", " min_x, min_y, min_z = 1000, 1000, 1000\n", " \n", " for c in r.coords:\n", " max_z = max(c[0], max_z)\n", " max_y = max(c[1], max_y)\n", " max_x = max(c[2], max_x)\n", " \n", " min_z = min(c[0], min_z)\n", " min_y = min(c[1], min_y)\n", " min_x = min(c[2], min_x)\n", " if (min_z == max_z or min_y == max_y or min_x == max_x or r.area > areas[-3]):\n", " for c in r.coords:\n", " segmented_ct_scan[c[0], c[1], c[2]] = 0\n", " else:\n", " index = (max((max_x - min_x), (max_y - min_y), (max_z - min_z))) / (min((max_x - min_x), (max_y - min_y) , (max_z - min_z)))" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "baa73057-81b1-803d-6cd6-5fbb1661df07" }, "source": [ "The `plot_3d` function plots the 3D numpy array of CT Scans. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "a639882e-626d-75a6-d3b6-103fc06e6f9f" }, "outputs": [], "source": [ "def plot_3d(image, threshold=-300):\n", " \n", " # Position the scan upright, \n", " # so the head of the patient would be at the top facing the camera\n", " p = image.transpose(2,1,0)\n", " p = p[:,:,::-1]\n", " \n", " verts, faces = measure.marching_cubes(p, threshold)\n", "\n", " fig = plt.figure(figsize=(10, 10))\n", " ax = fig.add_subplot(111, projection='3d')\n", "\n", " # Fancy indexing: `verts[faces]` to generate a collection of triangles\n", " mesh = Poly3DCollection(verts[faces], alpha=0.1)\n", " face_color = [0.5, 0.5, 1]\n", " mesh.set_facecolor(face_color)\n", " ax.add_collection3d(mesh)\n", "\n", " ax.set_xlim(0, p.shape[0])\n", " ax.set_ylim(0, p.shape[1])\n", " ax.set_zlim(0, p.shape[2])\n", "\n", " plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "a56d281a-816e-404d-c762-255b6cffa227" }, "outputs": [], "source": [ "plot_3d(segmented_ct_scan, 604)" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "c162afe0-b2b3-df59-d24d-c05900b39cf0" }, "source": [ "There are more methods that can be tried to get the regions of interest.\n", "\n", " 1. The left blood vessels can further be filtered using shape properties because we know that nodules are spherical in shape.\n", " 2. Super pixel segmentation can be further used and shape properties can be applied on segmented regions.\n", " 3. CNN architectures like UNet can also be used to generate candidate regions of interest. \n", "\n", "There are also a lot of good papers on Lung Segmentation and Nodule Candidate generation using Image Processing methods. I am posting the methods with a little summary of all of them:\n", "\n", " 1. [Automatic segmentation of lung nodules with growing neural gas and support vector machine][1]- The proposed method consists of the acquisition of computerized tomography images of the lung, the reduction of the volume of interest through techniques for the extraction of the thorax, extraction of the lung, and reconstruction of the original shape of the parenchyma. After that, growing neural gas (GNG) is applied to constrain even more the structures that are denser than the pulmonary parenchyma (nodules, blood vessels, bronchi, etc.). The next stage is the separation of the structures resembling lung nodules from other structures, such as vessels and bronchi. Finally, the structures are classified as either nodule or non-nodule, through shape and texture measurements together with support vector machine.\n", " 2. [Automated Segmentation of Lung Regions using Morphological Operators in CT scan][2]- The `get_segmented_lungs` method is a little modification of this paper. The threshold of 604 was taken from this paper.\n", " 3. [Pre-processing methods for nodule detection in lung CT][3]- This paper has used dot enhancement\n", "filter applied to the 3D matrix of voxel data. This 3D filter attempts to determine local geometrical characteristics for each voxel, computing the eigenvalues of the Hessian matrix and evaluating a \"likelihood\" function that was purposely built to discriminate between local morphology of linear, planar and spherical objects, modeled as having 3D Gaussian sections (Q. Li, S. Sone and K. Doi [6]). By applying this 3D filter to artificial images, we have verified the efficiency in detecting the Gaussian-like regions even in the cases were they are superimposed to non-Gaussian ones.\n", " 4. [Lung Nodule Detection using a Neural Classifier][4]: This paper discusses a dot-enhancement filter for nodule candidate selection and a neural classifier for false-positive finding reduction. The performance is evaluated as a fully automated computerized method for the detection of lung nodules in screening CT in the identification of lung cancers that may be missed during visual interpretation. \n", "\n", "Next I will be discussing about preprocessing of the LUNA16 dataset and using it to train the UNet model.\n", "\n", "\n", " [1]: http://www.sciencedirect.com/science/article/pii/S0010482512001412\n", " [2]: http://www.ijser.org/onlineResearchPaperViewer.aspx?Automated-Segmentation-of-Lung-Regions-using-Morphological-Operators-in-CT-scan.pdf\n", " [3]: https://arxiv.org/pdf/physics/0507153.pdf\n", " [4]: http://www.ijetch.org/papers/136-D058.pdf" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "a8dc2219-4d30-4dac-9460-80bab1216182" }, "source": [ "UNET for Candidate Point Generation\n", "-----------------------------------\n", "\n", "Nowadays, Deep Learning methods are achieving good results in Segmentation problems in medical imaging. One very famous architecture is UNET, which can be used for Nodule Candidates Points Generation in our case. Training of these networks is done using annotated datasets. The image processing methods used for candidate points generation above does not need any training data. We use the LUNA16 dataset for training our UNET model. \n", "\n", "In the LUNA16 dataset, each CT Scan is annotated with nodule points and the radius of the nodule which is used to generate the binary mask. I will first talk about preprocessing of the LUNA16 dataset. In the dataset, the CT Scans are saved in '.mhd' files and SimpleITK is used to read the image. I have defined three functions:\n", "\n", " - **load_itk** - Used to read a CT_Scan for the '.mhd' file.\n", " - **world_2_voxel**- Convert world coordinates to voxel coordinates.\n", " - **voxel_2_world**- Convert voxel coordinates to world coordinates." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "37acee08-d7f2-a666-2be4-b2e766cd6422" }, "outputs": [], "source": [ "'''\n", "This funciton reads a '.mhd' file using SimpleITK and return the image array, \n", "origin and spacing of the image.\n", "'''\n", "def load_itk(filename):\n", " # Reads the image using SimpleITK\n", " itkimage = sitk.ReadImage(filename)\n", " \n", " # Convert the image to a numpy array first and then shuffle the dimensions to get axis in the order z,y,x\n", " ct_scan = sitk.GetArrayFromImage(itkimage)\n", " \n", " # Read the origin of the ct_scan, will be used to convert the coordinates from world to voxel and vice versa.\n", " origin = np.array(list(reversed(itkimage.GetOrigin())))\n", " \n", " # Read the spacing along each dimension\n", " spacing = np.array(list(reversed(itkimage.GetSpacing())))\n", " \n", " return ct_scan, origin, spacing\n", "\n", "'''\n", "This function is used to convert the world coordinates to voxel coordinates using \n", "the origin and spacing of the ct_scan\n", "'''\n", "def world_2_voxel(world_coordinates, origin, spacing):\n", " stretched_voxel_coordinates = np.absolute(world_coordinates - origin)\n", " voxel_coordinates = stretched_voxel_coordinates / spacing\n", " return voxel_coordinates\n", "\n", "'''\n", "This function is used to convert the voxel coordinates to world coordinates using \n", "the origin and spacing of the ct_scan.\n", "'''\n", "def voxel_2_world(voxel_coordinates, origin, spacing):\n", " stretched_voxel_coordinates = voxel_coordinates * spacing\n", " world_coordinates = stretched_voxel_coordinates + origin\n", " return world_coordinates" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "a1d9f2ac-af87-cdec-64ce-48da9fc1e25a" }, "source": [ "After reading the 3D CT Scans, we will first segment the lungs and then generate the binary mask of nodule regions. This will be done by the `create_nodule_mask` function. The `draw_circle` function is used to mark the nodule regions in the binary mask. 'cands' are the list of nodule points with the radius given in the `annotation.csv` file of LUNA16 dataset. At the end we save the resized CT Scan with its segmented lungs and binary mask of nodules." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "686a2dd6-7767-51c5-f5ce-66c63ff8dcb9" }, "outputs": [], "source": [ "def seq(start, stop, step=1):\n", "\tn = int(round((stop - start)/float(step)))\n", "\tif n > 1:\n", "\t\treturn([start + step*i for i in range(n+1)])\n", "\telse:\n", "\t\treturn([])\n", "\n", "'''\n", "This function is used to create spherical regions in binary masks\n", "at the given locations and radius.\n", "'''\n", "def draw_circles(image,cands,origin,spacing):\n", "\t#make empty matrix, which will be filled with the mask\n", "\tRESIZE_SPACING = [1, 1, 1]\n", " image_mask = np.zeros(image.shape)\n", "\n", "\t#run over all the nodules in the lungs\n", "\tfor ca in cands.values:\n", "\t\t#get middel x-,y-, and z-worldcoordinate of the nodule\n", "\t\tradius = np.ceil(ca[4])/2\n", "\t\tcoord_x = ca[1]\n", "\t\tcoord_y = ca[2]\n", "\t\tcoord_z = ca[3]\n", "\t\timage_coord = np.array((coord_z,coord_y,coord_x))\n", "\n", "\t\t#determine voxel coordinate given the worldcoordinate\n", "\t\timage_coord = world_2_voxel(image_coord,origin,spacing)\n", "\n", "\t\t#determine the range of the nodule\n", "\t\tnoduleRange = seq(-radius, radius, RESIZE_SPACING[0])\n", "\n", "\t\t#create the mask\n", "\t\tfor x in noduleRange:\n", "\t\t\tfor y in noduleRange:\n", "\t\t\t\tfor z in noduleRange:\n", "\t\t\t\t\tcoords = world_2_voxel(np.array((coord_z+z,coord_y+y,coord_x+x)),origin,spacing)\n", "\t\t\t\t\tif (np.linalg.norm(image_coord-coords) * RESIZE_SPACING[0]) < radius:\n", "\t\t\t\t\t\timage_mask[np.round(coords[0]),np.round(coords[1]),np.round(coords[2])] = int(1)\n", "\t\n", "\treturn image_mask\n", "\n", "'''\n", "This function takes the path to a '.mhd' file as input and \n", "is used to create the nodule masks and segmented lungs after \n", "rescaling to 1mm size in all directions. It saved them in the .npz\n", "format. It also takes the list of nodule locations in that CT Scan as \n", "input.\n", "'''\n", "def create_nodule_mask(imagePath, maskPath, cands):\n", "\t#if os.path.isfile(imagePath.replace('original',SAVE_FOLDER_image)) == False:\n", "\timg, origin, spacing = load_itk(imagePath)\n", "\n", "\t#calculate resize factor\n", " RESIZE_SPACING = [1, 1, 1]\n", "\tresize_factor = spacing / RESIZE_SPACING\n", "\tnew_real_shape = img.shape * resize_factor\n", "\tnew_shape = np.round(new_real_shape)\n", "\treal_resize = new_shape / img.shape\n", "\tnew_spacing = spacing / real_resize\n", "\t\n", "\t#resize image\n", "\tlung_img = scipy.ndimage.interpolation.zoom(img, real_resize)\n", " \n", " # Segment the lung structure\n", "\tlung_img = lung_img + 1024\n", "\tlung_mask = segment_lung_from_ct_scan(lung_img)\n", "\tlung_img = lung_img - 1024\n", "\n", "\t#create nodule mask\n", "\tnodule_mask = draw_circles(lung_img,cands,origin,new_spacing)\n", "\n", "\tlung_img_512, lung_mask_512, nodule_mask_512 = np.zeros((lung_img.shape[0], 512, 512)), np.zeros((lung_mask.shape[0], 512, 512)), np.zeros((nodule_mask.shape[0], 512, 512))\n", "\n", "\toriginal_shape = lung_img.shape\t\n", "\tfor z in range(lung_img.shape[0]):\n", "\t\toffset = (512 - original_shape[1])\n", "\t\tupper_offset = np.round(offset/2)\n", "\t\tlower_offset = offset - upper_offset\n", "\n", "\t\tnew_origin = voxel_2_world([-upper_offset,-lower_offset,0],origin,new_spacing)\n", "\n", "\t\tlung_img_512[z, upper_offset:-lower_offset,upper_offset:-lower_offset] = lung_img[z,:,:]\n", "\t\tlung_mask_512[z, upper_offset:-lower_offset,upper_offset:-lower_offset] = lung_mask[z,:,:]\n", "\t\tnodule_mask_512[z, upper_offset:-lower_offset,upper_offset:-lower_offset] = nodule_mask[z,:,:]\n", "\n", " # save images. \n", "\tnp.save(imageName + '_lung_img.npz', lung_img_512)\n", "\tnp.save(imageName + '_lung_mask.npz', lung_mask_512)\n", "\tnp.save(imageName + '_nodule_mask.npz', nodule_mask_512)\n", " " ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "212658ac-21e5-51bd-a5da-bc5de1f92d3c" }, "source": [ "After preprocessing the dataset, the next thing is to train the model for segmentation. The model is written in keras in the `unet_model` function. It takes a 2D slice as input and returns a 2D slice of the same size as output. There are few things to be kept in mind while training\n", "\n", " - We wont use the slices that has no nodule region in the mask for training.\n", " - Dataset augmentation is very important because are nodules are generally circular or spherical in shape and are of different radius. \n", " - Since the nodule regions are very less, the dataset is skewed. Thus, we should weight the loss function accordingly.\n", " - The model may overfit on the training dataset. Thus, Dropout or Spatial Dropout are used to avoid overfitting." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "7e9734f2-9e28-a919-ef6c-e62fdd09aa00" }, "outputs": [], "source": [ "# change the loss function\n", "def dice_coef(y_true, y_pred):\n", " smooth = 1.\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", "\n", "\n", "def dice_coef_loss(y_true, y_pred):\n", "\treturn -dice_coef(y_true, y_pred)\n", "\n", "'''\n", "The UNET model is compiled in this function.\n", "'''\n", "def unet_model():\n", "\tinputs = Input((1, 512, 512))\n", "\tconv1 = Convolution2D(64, 3, 3, activation='relu', border_mode='same')(inputs)\n", "\tconv1 = Dropout(0.2)(conv1)\n", "\tconv1 = Convolution2D(64, 3, 3, activation='relu', border_mode='same')(conv1)\n", "\tpool1 = MaxPooling2D(pool_size=(2, 2))(conv1)\n", "\n", "\tconv2 = Convolution2D(128, 3, 3, activation='relu', border_mode='same')(pool1)\n", "\tconv2 = Dropout(0.2)(conv2)\n", "\tconv2 = Convolution2D(128, 3, 3, activation='relu', border_mode='same')(conv2)\n", "\tpool2 = MaxPooling2D(pool_size=(2, 2))(conv2)\n", "\n", "\tconv3 = Convolution2D(256, 3, 3, activation='relu', border_mode='same')(pool2)\n", "\tconv3 = Dropout(0.2)(conv3)\n", "\tconv3 = Convolution2D(256, 3, 3, activation='relu', border_mode='same')(conv3)\n", "\tpool3 = MaxPooling2D(pool_size=(2, 2))(conv3)\n", "\n", "\tconv4 = Convolution2D(512, 3, 3, activation='relu', border_mode='same')(pool3)\n", "\tconv4 = Dropout(0.2)(conv4)\n", "\tconv4 = Convolution2D(512, 3, 3, activation='relu', border_mode='same')(conv4)\n", "\tpool4 = MaxPooling2D(pool_size=(2, 2))(conv4)\n", "\n", "\tconv5 = Convolution2D(1024, 3, 3, activation='relu', border_mode='same')(pool4)\n", "\tconv5 = Dropout(0.2)(conv5)\n", "\tconv5 = Convolution2D(1024, 3, 3, activation='relu', border_mode='same')(conv5)\n", "\n", "\tup6 = merge([UpSampling2D(size=(2, 2))(conv5), conv4], mode='concat', concat_axis=1)\n", "\tconv6 = Convolution2D(512, 3, 3, activation='relu', border_mode='same')(up6)\n", "\tconv6 = Dropout(0.2)(conv6)\n", "\tconv6 = Convolution2D(512, 3, 3, activation='relu', border_mode='same')(conv6)\n", "\n", "\tup7 = merge([UpSampling2D(size=(2, 2))(conv6), conv3], mode='concat', concat_axis=1)\n", "\tconv7 = Convolution2D(256, 3, 3, activation='relu', border_mode='same')(up7)\n", "\tconv7 = Dropout(0.2)(conv7)\n", "\tconv7 = Convolution2D(256, 3, 3, activation='relu', border_mode='same')(conv7)\n", "\n", "\tup8 = merge([UpSampling2D(size=(2, 2))(conv7), conv2], mode='concat', concat_axis=1)\n", "\tconv8 = Convolution2D(128, 3, 3, activation='relu', border_mode='same')(up8)\n", "\tconv8 = Dropout(0.2)(conv8)\n", "\tconv8 = Convolution2D(128, 3, 3, activation='relu', border_mode='same')(conv8)\n", "\n", "\tup9 = merge([UpSampling2D(size=(2, 2))(conv8), conv1], mode='concat', concat_axis=1)\n", "\tconv9 = Convolution2D(64, 3, 3, activation='relu', border_mode='same')(up9)\n", "\tconv9 = Dropout(0.2)(conv9)\n", "\tconv9 = Convolution2D(64, 3, 3, activation='relu', border_mode='same')(conv9)\n", "\n", "\tconv10 = Convolution2D(1, 1, 1, activation='sigmoid')(conv9)\n", "\n", "\tmodel = Model(input=inputs, output=conv10)\n", "\tmodel.summary()\n", "\tmodel.compile(optimizer=Adam(lr=1e-5), loss=dice_coef_loss, metrics=[dice_coef])\n", "\n", "\treturn model" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "09429b58-990b-8d50-cc13-3211d5bedf26" }, "source": [ "## Classification of Lung Nodules ##\n", "\n", "The candidate regions generated still has a lot of noise. Thus we need a classifier to classify the candidates as either Nodule or Non-Nodule. This is also referred to as False Positive Reduction step. I am sharing a list of useful resources with a summary of all of them:\n", "\n", " 1. [Pulmonary Nodule Classification with Convolutional Neural Networks][1]: They use a CNN to predict\n", "whether the image contains a pulmonary lesion. As baselines, they also look at using SVM, kNN, and logistic regression to perform the same task. LIDC dataset is used for training and the classification results are very good.\n", " 2. [Pulmonary Nodule Classification with Deep Convolutional Neural Networks on Computed Tomography Images][2]: They design a deep convolutional neural networks method for nodule classification, which has an advantage of autolearning representation and strong generalization ability. A specified network structure for nodule images is proposed to solve the recognition of three types of nodules, that is, solid, semisolid, and ground glass opacity (GGO). Deep convolutional neural networks are trained by 62,492 regions-of-interest (ROIs) samples including 40,772 nodules and 21,720 nonnodules from the Lung Image Database Consortium (LIDC) database. Experimental results demonstrate the effectiveness of the proposed method in terms of sensitivity and overall accuracy and that it consistently outperforms the competing methods.\n", " 3. [Deep Learning for the Classification of Lung Nodules][3]: They have trained a CNN and tested the results on different datasets. In order to understand the features of lung nodules, they further construct new datasets, based on the combination of artificial geometric nodules and some transformations of the original images, as well as a stochastic nodule shape model.\n", " 4. [Lung Nodule Classification Based on Deep Convolutional Neural Networks][4]: In this work, they present a method for classifying lung nodules based on CNNs. Training is performed by balancing the mini-batches on each stochastic gradient descent (SGD) iteration to address the lack of nodule samples compared to background samples. They show that the method outperforms a base feature-engineering method using the same techniques for other stages of lung nodule detection, and show that CNNs obtain competitive results.\n", "\n", "\n", " [1]: http://cs231n.stanford.edu/reports2016/324_Report.pdf\n", " [2]: https://www.hindawi.com/journals/cmmm/2016/6215085/\n", " [3]: https://arxiv.org/abs/1611.06651\n", " [4]: http://www.liv.ic.unicamp.br/~juliomb/resources/2016-lnc-ciarp.pdf" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "767f1baa-b704-148c-5eb5-c8416f3fa01f" }, "source": [ "Due to computational limitations the size of each cannot very large too. I will be considering cubic voxels of size 36. In the dataset, we have 1187 nodule points, which is a very less number for training a deep network. Thus we can cut many voxels around a nodule center and increase the size of the dataset for training. We will also sample equal number of negative examples from the image for training. The code for preprocessing the LUNA16 dataset is below:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "65a1aa80-6f32-5086-5b4b-9cbee4411df9" }, "outputs": [], "source": [ "'''\n", "This function takes the numpy array of CT_Scan and a list of coords from\n", "which voxels are to be cut. The list of voxel arrays are returned. We keep the \n", "voxels cubic because per pixel distance is same in all directions.\n", "'''\n", "def get_patch_from_list(lung_img, coords, window_size = 10):\n", "\tshape = lung_img.shape\n", "\toutput = []\n", "\tlung_img = lung_img + 1024\n", "\tfor i in range(len(coords)):\n", "\t\tpatch = lung_img[coords[i][0] - 18: coords[i][0] + 18,\n", "\t\t\t\t\t\t coords[i][1] - 18: coords[i][1] + 18,\n", "\t\t\t\t\t\t coords[i][2] - 18: coords[i][2] + 18]\t\t\t \n", "\t\toutput.append(patch)\n", "\t\n", "\treturn output\n", "\n", "'''\n", "Sample a random point from the image and return the coordinates. \n", "'''\n", "def get_point(shape):\n", "\tx = randint(50, shape[2] - 50)\n", "\ty = randint(50, shape[1] - 50)\n", "\tz = randint(20, shape[0] - 20)\n", "\treturn np.asarray([z, y, x])\n", "\n", "'''\n", "This function reads the training csv file which contains the CT Scan names and\n", "location of each nodule. It cuts many voxels around a nodule and takes random points as \n", "negative samples. The voxels are dumped using pickle. It is to be noted that the voxels are \n", "cut from original Scans and not the masked CT Scans generated while creating candidate\n", "regions.\n", "'''\n", "def create_data(path, train_csv_path):\n", "\tcoords, trainY = [], []\n", "\t\n", "\twith open(train_csv_path, 'rb') as f:\n", "\t\tlines = f.readlines()\n", "\t\t\n", "\t\tfor line in lines:\n", "\t\t\trow = line.split(',')\n", "\t\t\t\n", "\t\t\tif os.path.isfile(path + row[0] + '.mhd') == False:\n", "\t\t\t\tcontinue\n", "\n", "\t\t\tlung_img = sitk.GetArrayFromImage(sitk.ReadImage(path + row[0] + '.mhd'))\n", "\n", "\t\t\tfor i in range(-5, 5, 3):\n", "\t\t\t\tfor j in range(-5, 5, 3):\n", "\t\t\t\t\tfor k in range(-2, 3, 2):\n", "\t\t\t\t\t\tcoords.append([int(row[3]) + k, int(row[2]) + j, int(row[1]) + i])\n", "\t\t\t\t\t\ttrainY.append(True)\n", "\n", "\t\t\tfor i in range(60):\t\t\t\n", "\t\t\t\tcoords.append(get_point(lung_img.shape))\n", "\t\t\t\ttrainY.append(False)\n", "\n", "\t\t\ttrainX = get_patch_from_list(lung_img, coords)\n", "\n", "\t\t\tprint (trainX[0].shape, len(trainX))\n", "\n", "\t\t\tpickle.dump(np.asarray(trainX), open('traindata_' + str(counter) + '_Xtrain.p', 'wb'))\n", "\t\t\tpickle.dump(np.asarray(trainY, dtype = bool), open('traindata_' + str(counter) + '_Ytrain.p', 'wb'))\n", "\n", "\t\t\tcounter = counter + 1\n", "\t\t\tcoords, trainY = [], []" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "248c2aa6-a0cb-1bd6-85e4-0b1e6ea6cc2f" }, "source": [ "In the processing, there are few things that must be tried:\n", "\n", " 1. Do not sample negative examples randomly from CT Scan. Select using the candiate points or the False Positives Generated, because it will train the model to discriminate non-nodule regions better.\n", " 2. Dataset Augmentation should be done while training because the nodules are symmetric regions and of varying sizes, thus dataset augmentation is very likely to help in such cases.\n", "\n", "Next is the script for training a classifier with 3D Convolutional Neural Networks. In the preprocessed pickle data files, randomly select 20% files and put it in a seperate folder as test set. The model is trained in keras. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "b7c3eb38-8858-a114-9565-c5fdc95a94a0" }, "outputs": [], "source": [ "from keras.models import Sequential\n", "from keras.layers import Dense, Dropout, Activation, Flatten\n", "from keras.layers import Convolution3D, MaxPooling3D\n", "from keras.utils import np_utils\n", "from keras import backend as K\n", "\n", "'''\n", "Creates a keras model with 3D CNNs and returns the model.\n", "'''\n", "def classifier(input_shape, kernel_size, pool_size):\n", "\tmodel = Sequential()\n", "\n", "\tmodel.add(Convolution3D(16, kernel_size[0], kernel_size[1], kernel_size[2],\n", "\t border_mode='valid',\n", "\t input_shape=input_shape))\n", "\tmodel.add(Activation('relu'))\n", "\tmodel.add(MaxPooling3D(pool_size=pool_size))\n", "\tmodel.add(Convolution2D(32, kernel_size[0], kernel_size[1], kernel_size[2]))\n", "\tmodel.add(Activation('relu'))\n", "\tmodel.add(MaxPooling3D(pool_size=pool_size))\n", "\tmodel.add(Convolution3D(64, kernel_size[0], kernel_size[1], kernel_size[2]))\n", "\tmodel.add(Activation('relu'))\n", "\tmodel.add(MaxPooling3D(pool_size=pool_size))\n", "\tmodel.add(Dropout(0.25))\n", "\n", "\tmodel.add(Flatten())\n", "\tmodel.add(Dense(512))\n", "\tmodel.add(Activation('relu'))\n", "\tmodel.add(Dropout(0.5))\n", "\tmodel.add(Dense(128))\n", "\tmodel.add(Activation('relu'))\n", "\tmodel.add(Dropout(0.5))\n", "\tmodel.add(Dense(2))\n", "\tmodel.add(Activation('softmax'))\n", "\n", "\treturn model\n", "\n", "def train_classifier(input_shape):\n", "\tmodel = classifier(input_shape, (3, 3, 3), (2, 2, 2))\n", "\tmodel.compile(loss='categorical_crossentropy',\n", "\t\t optimizer='adadelta',\n", "\t\t metrics=['accuracy'])\n", "\t'''\n", "\tRead the preprocessed datafiles chunk by chunnk and train the model on that batch (trainX, trainY) using:'''\n", "\tmodel.train_on_batch(trainX, trainY, sample_weight=None)\n", "\t'''The model can be trained on many epochs using for loops'''\n", "\t\n", "\t'''\n", "\tAFter training the dataset we test our model of the test dataset, read the test file chunk by chunk and \n", "\ttest on each chunk (trainX, trainY) using:'''\n", "\tprint (model.test_on_batch(trainX, trainY, sample_weight=None))\n", "\t'''The average of all of this can be taken to obtain the final test score.'''\n", "\t\n", "\t'''After testing save the model using'''\n", "\tmodel.save('my_model.h5')\n", " " ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "1990786a-009e-ca2a-ebd9-a54b42d6437f" }, "source": [ "This is a simple 3D CNN architecture for classification. There are few things that anyone should try to improve the accuracy:\n", "1) Use of multi-scale CNNs to capture more features for classification.\n", "2) Extract the features of fully connected layers and use xgboost on it to improve classification accuracy.\n", "\n", "In summary, we first generate candidates using Image Processing or UNet and then cut voxels around all candidate points and classify them. The points are clustered using DBScan algorithm and probability of a cluster is the mean of probabilities of all points in that cluster. The activation of the fully connected network of the classifier is extracted and this features is created for each CT Scan in our dataset. There is one more xgboost or SVM classifier on our dataset to classify the samples. We take the mean of the probabilities of this classifier as the final output probability.\n", "\n", "I am currently working on training a classifier on the LUNA16 dataset. I will make the results available with dataset augmentation and about making a sample submission on the kernel as soon as it I complete it." ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "2d879e9c-e6fc-e584-7e6c-e4c7f23a9475" }, "source": [ "I will post the dataset augmentation functions and Spatial Dropout very soon. I will post more methods as soon I implement them.\n", "\n", "After getting the regions of interest, 3D voxels can be cut around regions of interest and used for classification.\n", "\n", "**Please upvote or leave a comment, if you liked the tutorial.**" ] } ], "metadata": { "_change_revision": 0, "_is_fork": false, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.0" } }, "nbformat": 4, "nbformat_minor": 0 }