{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "4fae53f0-743a-e0e8-2638-d3a56d09eb70"
      },
      "source": [
        "# Data Science Bowl 2017\n",
        "\n",
        "The yearly tradition continues! - this time with another computer vision problem. we have to classify whether someone will be diagnosed with **lung cancer** at some point during the next year. We are given DICOM files, which is a format that is often used for medical scans. Using CT scans from 1400 patients in the training set, we have to build a model which can predict on the patients in the test set.\n",
        "\n",
        "**Shameless plug:** If you have any questions or want to discuss competitions/hardware/games/anything with other Kagglers, then join the [KaggleNoobs Slack channel](https://go o.gl/gGWFXe)! I feel like it could use a lot more users :)\n",
        "\n",
        "**And as always, if this helped you, some upvotes would be very much appreciated! :D**"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "814fd7c2-9570-f2cb-6187-fd396eb9891c"
      },
      "outputs": [],
      "source": [
        "import numpy as np # linear algebra\n",
        "import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n",
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "import os\n",
        "import glob\n",
        "%matplotlib inline\n",
        "p = sns.color_palette()\n",
        "\n",
        "os.listdir('../input')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "79562094-c8a2-a737-452a-4fc256563a6e"
      },
      "source": [
        "It looks like this time, we have been given the sample submission, the labels for the training set, and a sample submission for a test set.\n",
        "\n",
        "# The image data\n",
        "\n",
        "First, let's look at the sample images!"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "fe22bf62-1b6e-afd1-f7a4-22a52fc64904"
      },
      "outputs": [],
      "source": [
        "for d in os.listdir('../input/sample_images'):\n",
        "    print(\"Patient '{}' has {} scans\".format(d, len(os.listdir('../input/sample_images/' + d))))\n",
        "print('----')\n",
        "print('Total patients {} Total DCM files {}'.format(len(os.listdir('../input/sample_images')), \n",
        "                                                      len(glob.glob('../input/sample_images/*/*.dcm'))))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "c40cf73c-3678-4bfd-2984-64bc7c0b69ef"
      },
      "outputs": [],
      "source": [
        "patient_sizes = [len(os.listdir('../input/sample_images/' + d)) for d in os.listdir('../input/sample_images')]\n",
        "plt.hist(patient_sizes, color=p[2])\n",
        "plt.ylabel('Number of patients')\n",
        "plt.xlabel('DICOM files')\n",
        "plt.title('Histogram of DICOM count per patient')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "dcd6ddf2-9d4b-cd55-0c84-dd1003530860"
      },
      "source": [
        "We can see that the sample_images directory is made up of a bunch of subdirectories, each representing a single patient ID and containing about 100-300 DICOM files inside (except one with over 400).\n",
        "\n",
        "What about file size?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "06820c48-af7f-ea3d-8c69-b780b1952648"
      },
      "outputs": [],
      "source": [
        "sizes = [os.path.getsize(dcm)/1000000 for dcm in glob.glob('../input/sample_images/*/*.dcm')]\n",
        "print('DCM file sizes: min {:.3}MB max {:.3}MB avg {:.3}MB std {:.3}MB'.format(np.min(sizes), \n",
        "                                                       np.max(sizes), np.mean(sizes), np.std(sizes)))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "5beb4fab-3df1-6487-ffef-cbdbb0b80cd2"
      },
      "source": [
        "We can see here that basically all the DICOM files are the same 0.525MB size - which makes it seem like they are uncompressed files with the same amount of information.\n",
        "\n",
        "We'll come back to images to load them in a second, but first let's look at the other two files.\n",
        "\n",
        "# Training set"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "4c2798d2-d22e-2ae8-a1a4-17b969d0be16"
      },
      "outputs": [],
      "source": [
        "df_train = pd.read_csv('../input/stage1_labels.csv')\n",
        "df_train.head()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "e1ea9231-52d4-a914-e12f-95508365051a"
      },
      "source": [
        "Simple enough. Let's look at the data stats:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "0bd13534-5177-630e-8ff5-8c80219f13d4"
      },
      "outputs": [],
      "source": [
        "print('Number of training patients: {}'.format(len(df_train)))\n",
        "print('Cancer rate: {:.4}%'.format(df_train.cancer.mean()*100))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "f1b993aa-0815-7ff7-640a-c8c1f7a4820d"
      },
      "source": [
        "Here we can see that in the training set we have a total of just under 1400 patients, of which just over a quarter have had lung cancer diagnosed in the year after these scans were taken.\n",
        "\n",
        "Let's see if we can exploit this information to beat the benchmark on the leaderboard!\n",
        "\n",
        "### Naive Submission\n",
        "\n",
        "Since the evaluation metric used in this competition is [LogLoss](https://www.kaggle.com/wiki/LogarithmicLoss) and not something like [AUC](https://www.kaggle.com/wiki/AreaUnderCurve), this means that we can often gain an improvement just by aligning the probabilities of our sample submission to that of the training set.\n",
        "\n",
        "Before I try making a naive submission, I will calculate what the score of this submission would be on the training set to get a comparison."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "24e3fbbc-f244-4837-801d-f0336a9f4282"
      },
      "outputs": [],
      "source": [
        "from sklearn.metrics import log_loss\n",
        "logloss = log_loss(df_train.cancer, np.zeros_like(df_train.cancer) + df_train.cancer.mean())\n",
        "print('Training logloss is {}'.format(logloss))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "d964907d-989f-9f3e-176f-63e5c0404762"
      },
      "source": [
        "Wow, not bad. Let's turn it into a submission."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "2ec6faae-4498-8dd3-747f-166b852eb08d"
      },
      "outputs": [],
      "source": [
        "sample = pd.read_csv('../input/stage1_sample_submission.csv')\n",
        "sample['cancer'] = df_train.cancer.mean()\n",
        "sample.to_csv('naive_submission.csv', index=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "7b4e658b-8f6b-6df5-10d2-d219a2d0586b"
      },
      "source": [
        "This submission scores **0.60235** on the leaderboard - you can try it out for yourself by heading over to the **Output** tab at the top of the kernel.\n",
        "\n",
        "The fact that the score is worse here shows us that we have [overfitted](https://xkcd.com/605/) to our training data. The mean of the test set is different than the mean of the training set, and while this may only be a small difference, it is the reason why the score is worse on the leaderboard. But we won't be winning any prizes for this submission anyway, so it's time to move on :)\n",
        "\n",
        "### Leakbusters\n",
        "\n",
        "Just as a quick sanity check, I'm going to check to see whether there's any observable relationship between Patient ID and whether they have cancer or not. Let's hope they used better random seeds than TalkingData #neverforget\n",
        "\n",
        "Note that in the following graph, it is also in order of PatientID as the rows have been presorted."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "9dca4e7f-b3e2-7dd3-5563-12c611b18f8e"
      },
      "outputs": [],
      "source": [
        "targets = df_train['cancer']\n",
        "plt.plot(pd.rolling_mean(targets, window=10), label='Sliding Window 10')\n",
        "plt.plot(pd.rolling_mean(targets, window=50), label='Sliding Window 50')\n",
        "plt.xlabel('rowID')\n",
        "plt.ylabel('Mean cancer')\n",
        "plt.title('Mean target over rowID - sliding mean')\n",
        "plt.legend()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "51f6ddcc-8a6f-aaac-ff20-72f166837f84"
      },
      "outputs": [],
      "source": [
        "print('Accuracy predicting no cancer: {}%'.format((df_train['cancer'] == 0).mean()))\n",
        "print('Accuracy predicting with last value: {}%'.format((df_train['cancer'] == df_train['cancer'].shift()).mean()))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "9fefc4c0-a0d6-152b-1be6-0f7676ad4efb"
      },
      "source": [
        "There no leak here which is immediately apparent to me - it looks well sorted. Feel free to disagree if you've noticed something!"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "97650778-09fb-bc59-607d-bdabae2eb749"
      },
      "source": [
        "# Test set (Stage 1)\n",
        "\n",
        "After looking at the training file for a bit, let's take a brief look at the test set/sample submission."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "4c67128f-d764-82b4-14f3-9035b50bb235"
      },
      "outputs": [],
      "source": [
        "sample = pd.read_csv('../input/stage1_sample_submission.csv')\n",
        "sample.head()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "1f2c3c95-cda7-107f-3daf-57f3b57beaff"
      },
      "outputs": [],
      "source": [
        "print('The test file has {} patients'.format(len(sample)))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "e1dd4a15-97cc-cb7e-5ce3-79be437604df"
      },
      "source": [
        "Nothing out of the ordinary here, the submission file is arranged very similarly to the training csv, except that we submit a probability instead of a class predction. There are less than 200 samples, which means that we need to watch out for overfitting on the leaderboard.\n",
        "\n",
        "It is actually possible to get a perfect score only through brute force in Stage 1 - as we have 200 samples and over 200 submissions to test on the LB. So expect lots of overfitting. [Trust no one, not even yourself.](http://i1.kym-cdn.com/entries/icons/original/000/017/046/BptVE1JIEAAA3dT.jpg)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "e34fbeb1-fc2b-01cf-9067-4ca44f254490"
      },
      "source": [
        "# DICOMs\n",
        "\n",
        "The dicom package was finally added! Here's some quick exploration of the files and what they contain.\n",
        "\n",
        "First, I'll just load a random image from the lot."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "3b9cda40-c2b5-35f9-7e32-6f776202b32a"
      },
      "outputs": [],
      "source": [
        "import dicom"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "d7fc983f-0eca-f256-2c4d-3e06048a68a8"
      },
      "outputs": [],
      "source": [
        "dcm = '../input/sample_images/0a38e7597ca26f9374f8ea2770ba870d/4ec5ef19b52ec06a819181e404d37038.dcm'\n",
        "print('Filename: {}'.format(dcm))\n",
        "dcm = dicom.read_file(dcm)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "4ca1a85d-659e-fa46-25a2-2c9b7d2cb271"
      },
      "source": [
        "Now I've read a random image, let's view the information inside."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "d6752daa-099a-fbce-f7fe-4ce894aa5d88"
      },
      "outputs": [],
      "source": [
        "dcm"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "61fd340c-bfe9-f30c-a4fe-8d5db93d7ea5"
      },
      "source": [
        "It looks this data has been intentionally anonymised in order to keep this a computer vision problem. Notably, the birth date has been anonymised to January 1st, 1900. Age could otherwise be an important feature for predicting lung cancer.\n",
        "\n",
        "There are two things here that I think are significant, slice location (this sounds like it could be the z-position of the scan?) and the 'Pixel Data'.\n",
        "\n",
        "We can retrieve a image as a **numpy array** by calling `dcm.pixel_array`, and we can then replace the -2000s, which are essentially NAs, with 0s (thanks r4m0n)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "f6ad63f4-7443-4107-e3bd-9d19d53bab3e"
      },
      "outputs": [],
      "source": [
        "img = dcm.pixel_array\n",
        "img[img == -2000] = 0\n",
        "\n",
        "plt.axis('off')\n",
        "plt.imshow(img)\n",
        "plt.show()\n",
        "\n",
        "plt.axis('off')\n",
        "plt.imshow(-img) # Invert colors with -\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "73196413-2230-505a-614c-f44e9dc7aa75"
      },
      "source": [
        "Let's plot a few more images at random:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "aa949b02-1267-eb26-3877-1a9dcd849cbb"
      },
      "outputs": [],
      "source": [
        "def dicom_to_image(filename):\n",
        "    dcm = dicom.read_file(filename)\n",
        "    img = dcm.pixel_array\n",
        "    img[img == -2000] = 0\n",
        "    return img"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "9815c2df-b19d-8d10-9d12-d5c46c388aa4"
      },
      "outputs": [],
      "source": [
        "files = glob.glob('../input/sample_images/*/*.dcm')\n",
        "\n",
        "f, plots = plt.subplots(4, 5, sharex='col', sharey='row', figsize=(10, 8))\n",
        "for i in range(20):\n",
        "    plots[i // 5, i % 5].axis('off')\n",
        "    plots[i // 5, i % 5].imshow(dicom_to_image(np.random.choice(files)), cmap=plt.cm.bone)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "e9442584-8cf2-f1ed-3f19-f9743feb4220"
      },
      "source": [
        "This gives us some idea with the sort of images we're dealing with. Now, let's try to reconstruct the layers of the body from which the images were taken, by taking a single patient and sorting his scans by Slice Location."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "e63b4765-28cb-ae97-c638-819e9fb3b038"
      },
      "outputs": [],
      "source": [
        "def get_slice_location(dcm):\n",
        "    return float(dcm[0x0020, 0x1041].value)\n",
        "\n",
        "# Returns a list of images for that patient_id, in ascending order of Slice Location\n",
        "def load_patient(patient_id):\n",
        "    files = glob.glob('../input/sample_images/{}/*.dcm'.format(patient_id))\n",
        "    imgs = {}\n",
        "    for f in files:\n",
        "        dcm = dicom.read_file(f)\n",
        "        img = dcm.pixel_array\n",
        "        img[img == -2000] = 0\n",
        "        sl = get_slice_location(dcm)\n",
        "        imgs[sl] = img\n",
        "        \n",
        "    # Not a very elegant way to do this\n",
        "    sorted_imgs = [x[1] for x in sorted(imgs.items(), key=lambda x: x[0])]\n",
        "    return sorted_imgs"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "c0309cb1-e9c3-efd2-e859-b5e4d7a096af"
      },
      "outputs": [],
      "source": [
        "pat = load_patient('0a38e7597ca26f9374f8ea2770ba870d')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "6a1db123-a677-477a-8428-6c1eae44c62a"
      },
      "source": [
        "Now that we have the images of a patient sorted by position in the body, we can plot them to see how this varies.\n",
        "\n",
        "It's worth noting that this patient does **not** have cancer.\n",
        "\n",
        "### Sorted Slices of Patient 0a38e7597ca26f9374f8ea2770ba870d - No cancer"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "c0e0d494-36e2-73e3-63e9-dffb88a27b9d"
      },
      "outputs": [],
      "source": [
        "f, plots = plt.subplots(11, 10, sharex='all', sharey='all', figsize=(10, 11))\n",
        "# matplotlib is drunk\n",
        "#plt.title('Sorted Slices of Patient 0a38e7597ca26f9374f8ea2770ba870d - No cancer')\n",
        "for i in range(110):\n",
        "    plots[i // 10, i % 10].axis('off')\n",
        "    plots[i // 10, i % 10].imshow(pat[i], cmap=plt.cm.bone)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "210fcd93-ba72-ef3c-8027-b17c31762746"
      },
      "source": [
        "### Sorted Slices of Patient 0acbebb8d463b4b9ca88cf38431aac69 - Cancer"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "ef841809-280c-5c9b-c009-dcaf1afcd325"
      },
      "outputs": [],
      "source": [
        "pat = load_patient('0acbebb8d463b4b9ca88cf38431aac69')\n",
        "f, plots = plt.subplots(21, 10, sharex='all', sharey='all', figsize=(10, 21))\n",
        "for i in range(203):\n",
        "    plots[i // 10, i % 10].axis('off')\n",
        "    plots[i // 10, i % 10].imshow(pat[i], cmap=plt.cm.bone)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "d6aae6a8-c56c-900c-7023-c2b444bbb936"
      },
      "source": [
        "Okay, so it looks like my theory that \"Slice Position\" refers to the z-position of the scan was correct (it probably said this in the documentation but who reads documentation??).\n",
        "\n",
        "We can actually use this to reconstruct a 3D model of the of the torso by simply concatenating the images together. Then something like a 3D convolutional network could be applied on top in order to identify points of interest in 3D space. Interesting stuff.\n",
        "\n",
        "I'm going to try to make an animated gif/video now, which should show this better than a bunch of plots. I've so far found no way to do this in python :("
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "c087bd50-2380-30d6-d571-22591f263447"
      },
      "outputs": [],
      "source": [
        "# This function takes in a single frame from the DICOM and returns a single frame in RGB format.\n",
        "def normalise(img):\n",
        "    normed = (img / 14).astype(np.uint8) # Magic number, scaling to create int between 0 and 255\n",
        "    img2 = np.zeros([*img.shape, 3], dtype=np.uint8)\n",
        "    for i in range(3):\n",
        "        img2[:, :, i] = normed\n",
        "    return img2"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "8dad2388-4e2c-6876-a189-22ab0fc99ba8"
      },
      "outputs": [],
      "source": [
        "npat = [normalise(p) for p in pat]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "6d25c206-03b1-ac6b-51d8-3a26c037bbfa"
      },
      "outputs": [],
      "source": [
        "pat = load_patient('0acbebb8d463b4b9ca88cf38431aac69')\n",
        "\n",
        "import matplotlib.animation as animation\n",
        "def animate(pat, gifname):\n",
        "    # Based on @Zombie's code\n",
        "    fig = plt.figure()\n",
        "    anim = plt.imshow(pat[0], cmap=plt.cm.bone)\n",
        "    def update(i):\n",
        "        anim.set_array(pat[i])\n",
        "        return anim,\n",
        "    \n",
        "    a = animation.FuncAnimation(fig, update, frames=range(len(pat)), interval=50, blit=True)\n",
        "    a.save(gifname, writer='imagemagick')\n",
        "    \n",
        "animate(pat, 'test.gif')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "6f9ced35-1d10-62b0-b281-f7d3f42bdabb"
      },
      "outputs": [],
      "source": [
        "IMG_TAG = \"\"\"<img src=\"data:image/gif;base64,{0}\">\"\"\"\n",
        "\n",
        "import base64\n",
        "from IPython.display import HTML\n",
        "\n",
        "def display_gif(fname):\n",
        "    data = open(fname, \"rb\").read()\n",
        "    data = base64.b64encode(data)\n",
        "    return HTML(IMG_TAG.format(data))\n",
        "\n",
        "display_gif(\"test.gif\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "04a98a7d-b93b-dcd7-8f1d-e330bebe2e08"
      },
      "source": [
        "The gif doesn't show up in the notebook, but if you want to see it head over to the 'Output' tab and you can see the gif playing!"
      ]
    }
  ],
  "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
}