{"cells":[{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"import os\nimport random\n\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nimport matplotlib.pyplot as plt\n\nimport pydicom\n\nfrom skimage.measure import regionprops, label\nfrom skimage.segmentation import mark_boundaries\n\nfrom tqdm import tqdm_notebook","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"# Get directory names/locations\ndata_root = os.path.abspath(\"../input/rsna-intracranial-hemorrhage-detection/\")\n\ntrain_img_root = data_root + \"/stage_1_train_images/\"\ntest_img_root  = data_root + \"/stage_1_test_images/\"\n\ntrain_labels_path = data_root + \"/stage_1_train.csv\"\ntest_labels_path  = data_root + \"/stage_1_test.csv\"\n\n# Create list of paths to actual training data\ntrain_img_paths = os.listdir(train_img_root)\ntest_img_paths  = os.listdir(test_img_root)\n\n# Dataset size\nnum_train = len(train_img_paths)\nnum_test  = len(test_img_paths)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"def create_efficient_df(data_path):\n    \n    # Define the datatypes we're going to use\n    final_types = {\n        \"ID\": \"str\",\n        \"Label\": \"float16\"\n    }\n    features = list(final_types.keys())\n    \n    # Use chunks to import the data so that less efficient machines can only use a \n    # specific amount of chunks on import\n    df_list = []\n\n    chunksize = 1_000_000\n\n    for df_chunk in pd.read_csv(data_path, dtype=final_types, chunksize=chunksize): \n        df_list.append(df_chunk)\n        \n    df = pd.concat(df_list)\n    df = df[~df.isin([np.nan, np.inf, -np.inf]).any(1)]\n\n    del df_list\n\n    return df\n\ntrain_labels_df = create_efficient_df(train_labels_path)\ntrain_labels_df[train_labels_df[\"Label\"] > 0].head()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"hem_types = [\n    \"epidural\",\n    \"intraparenchymal\",\n    \"intraventricular\",\n    \"subarachnoid\",\n    \"subdural\",\n    \"any\"\n]\n\nnew_cols = [\n    \"id\",\n    \"type_0\",\n    \"type_1\",\n    \"type_2\",\n    \"type_3\",\n    \"type_4\",\n    \"type_5\"\n]\n\nnum_ids = int(train_labels_df.shape[0] / len(hem_types))\nprint(\"Number of unique patient IDs: {}\".format(num_ids))\n\nempty_array = np.ones((num_ids, len(new_cols)))\nhem_df = pd.DataFrame(data=empty_array, columns=new_cols)\n\n# Fill in the ID of each image\nhem_df[\"id\"] = list(train_labels_df.iloc[::len(hem_types)][\"ID\"].str.split(pat=\"_\").str[1])\n    \n# Fill in the categorical columns of each image\nfor hem_ix, hem_col in enumerate(list(hem_df)[1:]):\n    hem_df[hem_col] = list(train_labels_df.iloc[hem_ix::len(hem_types), 1])\n    \nhem_df.info()\nhem_df[hem_df[\"type_5\"] > 0].head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","collapsed":true,"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":false},"cell_type":"markdown","source":"# Introduction\n\nHello, my name is Mick. I've written a couple of other notebooks for this competition but today I'm going to actually dive into some image analysis. This notebook is meant to be a demonstration of some image analysis techniques and as such will mostly be operating on only one or two images, rather than creating a pipeline that analyzes many images in a rapid manner. The main techniques I'll be focusing on will revolve around the Sci-Kit Image library and use segmentation methods to split the images into their component parts. Our goal will be something like pulling out all of the regions that are likely to contain hemorrhages within an individual CT scan. \n\nIf you want an introduction to the [file format these images come in](https://www.kaggle.com/smit2300/dicom-patient-analysis), [some anatomical background on hemorrhages](https://www.kaggle.com/smit2300/hemorrhage-medical-introduction), or otherwise [information on exploring this dataset](https://www.kaggle.com/smit2300/hemorrhage-eda-encoding-dicom-introduction), you can check out any of my other notebooks in this competition to get a feel for the RSNA's Intracranial Hemorrhage dataset. \n\nNow that we're settled, lets start doing some image analysis!"},{"metadata":{},"cell_type":"markdown","source":"# Image Viewing\n\nThe first thing I'll do is simply pick out and image that we can segment and try to find some interesting features of the image."},{"metadata":{"trusted":true},"cell_type":"code","source":"CERTAINTY = 0.95\n\n# Filter the dataframe to search for epidural hemorrhages\nepi_df = hem_df[(hem_df[\"type_0\"] > CERTAINTY) & (hem_df[\"type_1\"] < CERTAINTY) & (hem_df[\"type_2\"] < CERTAINTY) & (hem_df[\"type_3\"] < CERTAINTY) & (hem_df[\"type_4\"] < CERTAINTY)]\n\n# Custom indices of images that contain good looking hemorrhages to me (please suggest better image if anyone know of any!)\nepi_ix = 6\n\n# Slice out the record at the chosen index\nepi_record = epi_df.iloc[epi_ix, :]\n\n# Get the image path from the record\nepi_path = train_img_root + \"ID_\" + epi_record[\"id\"] + \".dcm\"\n\n# Use PyDICOM to open the image and get array data\nepidural_frame = pydicom.dcmread(epi_path).pixel_array\n\n# Normalize the array between 0 and 255\nepidural_frame = np.interp(epidural_frame, (epidural_frame.min(), epidural_frame.max()), (0, 255))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_frame(img_array, cmap, title):\n    plt.figure(figsize=(8,8))\n    plt.imshow(img_array, cmap=cmap)\n    plt.title(title, fontsize=16)\n    plt.axis(\"off\")\n    \nplot_frame(epidural_frame, \"bone\", \"CT Scan of Epidural Hemorrhage\")\nplot_frame(epidural_frame, \"hot\", \"CT Scan of Epidural Hemorrhage\")\nplot_frame(epidural_frame, \"viridis\", \"CT Scan of Epidural Hemorrhage\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Statistical Brightness Analysis\n\nAlright we've found an example of each type of exclusive hemorrhage. There aren't too many explicit differences between the images and I'm actually not even particularly certain that the images I've chosen are good representations of each type of hemorrhage (If anyone has suggestions for good exemplary images then please feel free to leave a comment)! However, we can see some clear artifacts in these images that I'll show you how to separate out from the image.\n\nThe first thing I'll do is check out a histogram for the brightness levesl of each image. We may be able to find regions that are representative of each of the types of hemorrhage and brightness that might indicate hemorrhage type. For the next couple steps I'll show analysis of just an epidural image and then we can expand that to encompass the other hemorrhage types."},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(epidural_frame.flatten())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see some clear peaks on this histogram for image brightness. The regions I'm going to focus on are those of high density. Since we're going to mostly be looking for hemorrhages in this dataset we want to analyze the regions of the image that correspond to the skull and any perturbations it might be seeing."},{"metadata":{},"cell_type":"markdown","source":"# Image Masking\n\nNow that we know the distribution of brightness within the image we can create some masks that capture specific regions of this histogram. This is a fairly primitive method to segment an image and many algorithms exist for adaptively filtering the image, but for explanatory purposes I'll just be using a hard threshold of ~200 to segment these images.\n\nLet's view the image, the binary mask that is created from thresholding, and the result of multiplying the original image with the mask."},{"metadata":{"trusted":true},"cell_type":"code","source":"dense_mask   = (epidural_frame > 200).astype(int)\ndense_frame  = (epidural_frame * dense_mask).astype(int)\ndense_coords = np.argwhere(dense_frame)\n\nplot_frame(epidural_frame, \"viridis\", \"Original Epidural Hemorrhage Image\")\nplot_frame(dense_mask, \"hot\", \"Region of High Density Segmentation Mask\")\nplot_frame(dense_frame, \"viridis\", \"Regions of High Density\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"norm_frame = epidural_frame / 255.0\nsegmented_frame = mark_boundaries(norm_frame, dense_mask, color=(255,0,0))\n\nplot_frame(segmented_frame, \"bone\", \"Marked Boundaries of High Density\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Awesome! We can see our image clearly segmented to specifically highlight regions of the highest density on the CT scan. This methodology has been used in slightly augmented forms for years and is sometimes used to create datasets for deep learning. \n\nThe image that we've chosen seems to have a small artifact on what I believe is the anterior left side of the patient's skull. There is also a discontinuity in density that wasn't immediately clear in the original scan images. It's possible that we've found an anomalous region of the image that could be the result of injury and thus a hemorrhage. For later notebooks I'll try to analyze regions of specific injury and other stops along the way for deep learning models, but for this notebook the last thing I want to show is the Scikit-Image `regionprops()` method."},{"metadata":{},"cell_type":"markdown","source":"# regionprops()\n\nThe regionprops() method may be familiar if you come from a MATLAB background. The method employs an algorithm that finds all regions within an image that meet a set of labelled criteria provided by the programmer. We've created that criteria by making our boolean mask earlier. In order to pass the labelling image to the `regionprops()` method we first have to turn it into a Scikit-Image safe labelling array using the `label()` method. \n\nOnce we've used the `regionprops()` method to split the image into specific regions we can iterate over those and print out information on all of the continuous regions found within the image. This is a useful process for things like computing area of shapes or performing morphological calculations on regions within an image."},{"metadata":{"trusted":true},"cell_type":"code","source":"label_img = label(dense_mask)\nregions   = regionprops(label_image=label_img, intensity_image=dense_frame)\n\nlargest_area = max([regions[x].area for x in range(len(regions))])\nprint(\"Largest high density area: {}\".format(largest_area))\n\nnum_bones = len(regions)\nprint(\"Number of high density regions: {}\".format(num_bones))\n\nfor prop_ix, props in enumerate(regions):\n    print(\"\\nHotspot number {}:\".format(prop_ix+1))\n    print(\"Area: %d\" % (props.area))\n    print(\"Centroid: (%.2f, %.2f)\" % (props.centroid[0],props.centroid[1]))\n    print(\"Mean Intensity: %.2f\" % (props.mean_intensity))\n    print(\"Max Intensity: %d\" % (props.max_intensity))\n    ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Conclusion\n\nThis has been a basic introduction to get started with basic image segmentation using Scikit-Image. If there are any methods that people would like to see in this notebook please feel free to leave a comment and I would be happy to include them in the notebook or learn about them if they're not something I'm already familiar with.\n\nI hope you enjoyed and maybe even learned a little from this notebook! I'll see you in next one."},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":1}