{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"
\n\n
\n\n

Mayo Clinic - STRIP AI - EDA

\n
CREATED BY: DARIEN SCHETTLER
\n\n
\n\n---\n\n
\n\n
\n πŸ›‘   WARNING:

THIS IS A WORK IN PROGRESS
\n
\n\n\n
\n πŸ‘   IF YOU FORK THIS OR FIND THIS HELPFUL   πŸ‘

PLEASE UPVOTE!

This was a lot of work for me and while it may seem silly, it makes me feel appreciated when others like my work. πŸ˜…\n
\n\n\n","metadata":{}},{"cell_type":"markdown","source":"

\n\n

\n\n

TABLE OF CONTENTS

\n\n---\n\n

0    IMPORTS

\n\n---\n\n

1    BACKGROUND INFORMATION

\n\n---\n\n

2    SETUP

\n\n---\n\n

3    HELPER FUNCTIONS

\n\n---\n\n

4    DATASET EXPLORATION

\n\n---\n\n

5    MODELLING

\n\n---","metadata":{}},{"cell_type":"markdown","source":"
\n\n\n\n

0  IMPORTS    ⤒

","metadata":{}},{"cell_type":"code","source":"!sudo apt-get update\n!sudo apt-get -y install libvips-dev\n!pip install pyvips\nimport pyvips","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"\\n... IMPORTS STARTING ...\\n\")\n\nprint(\"\\n\\tVERSION INFORMATION\")\n# Machine Learning and Data Science Imports\nimport tensorflow as tf; print(f\"\\t\\t– TENSORFLOW VERSION: {tf.__version__}\");\nimport tensorflow_hub as tfhub; print(f\"\\t\\t– TENSORFLOW HUB VERSION: {tfhub.__version__}\");\nimport tensorflow_addons as tfa; print(f\"\\t\\t– TENSORFLOW ADDONS VERSION: {tfa.__version__}\");\nimport tensorflow_io as tfio; print(f\"\\t\\t– TENSORFLOW I/O VERSION: {tfio.__version__}\");\nimport pandas as pd; pd.options.mode.chained_assignment = None;\nimport numpy as np; print(f\"\\t\\t– NUMPY VERSION: {np.__version__}\");\nimport sklearn; print(f\"\\t\\t– SKLEARN VERSION: {sklearn.__version__}\");\nfrom sklearn.preprocessing import RobustScaler, PolynomialFeatures\nfrom pandarallel import pandarallel; pandarallel.initialize();\nfrom sklearn.model_selection import GroupKFold, StratifiedKFold\nfrom scipy.spatial import cKDTree\n\n# # RAPIDS\n# import cudf, cupy, cuml\n# from cuml.neighbors import NearestNeighbors\n# from cuml.manifold import TSNE, UMAP\n\n# Built In Imports\nfrom kaggle_datasets import KaggleDatasets\nfrom collections import Counter\nfrom datetime import datetime\nfrom glob import glob\nimport openslide\nimport warnings\nimport requests\nimport hashlib\nimport imageio\nimport IPython\nimport sklearn\nimport urllib\nimport zipfile\nimport pickle\nimport random\nimport shutil\nimport string\nimport json\nimport math\nimport time\nimport gzip\nimport ast\nimport sys\nimport io\nimport os\nimport gc\nimport re\n\n# Visualization Imports\nfrom matplotlib.colors import ListedColormap\nfrom matplotlib.patches import Rectangle\nimport matplotlib.patches as patches\nimport plotly.graph_objects as go\nimport matplotlib.pyplot as plt\nfrom tqdm.notebook import tqdm; tqdm.pandas();\nimport plotly.express as px\nimport tifffile as tif\nimport seaborn as sns\nfrom PIL import Image, ImageEnhance; Image.MAX_IMAGE_PIXELS = 5_000_000_000;\nimport matplotlib; print(f\"\\t\\t– MATPLOTLIB VERSION: {matplotlib.__version__}\");\nfrom matplotlib import animation, rc; rc('animation', html='jshtml')\nimport plotly\nimport PIL\nimport cv2\n\nimport plotly.io as pio\nprint(pio.renderers)\n\ndef seed_it_all(seed=7):\n \"\"\" Attempt to be Reproducible \"\"\"\n os.environ['PYTHONHASHSEED'] = str(seed)\n random.seed(seed)\n np.random.seed(seed)\n tf.random.set_seed(seed)\n\n \nprint(\"\\n\\n... IMPORTS COMPLETE ...\\n\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n\n\n\n

1  BACKGROUND INFORMATION    ⤒

\n\n---\n","metadata":{}},{"cell_type":"markdown","source":"
\n\n

1.1 BASIC COMPETITION INFORMATION

\n\n---\n\n
PRIMARY TASK DESCRIPTION\n\nThe goal of this competition is to classify the blood clot origins in ischemic stroke. Using whole slide digital pathology images, you'll build a model that differentiates between the two major acute ischemic stroke (AIS) etiology subtypes: cardiac and large artery atherosclerosis.\n\n
BASIC BACKGROUND INFORMATION\n\nStroke remains the second-leading cause of death worldwide. Each year in the United States, over 700,000 individuals experience an ischemic stroke caused by a blood clot blocking an artery to the brain. A second stroke (23% of total events are recurrent) worsens the chances of the patient’s survival. However, subsequent strokes may be mitigated if physicians can determine stroke etiology, which influences the therapeutic management following stroke events.\n\nDuring the last decade, mechanical thrombectomy has become the standard of care treatment for acute ischemic stroke from large vessel occlusion. \n\n
\n\n
\n\nAs a result, retrieved clots became amenable to analysis. \n* Healthcare professionals are currently attempting to apply deep learning-based methods to predict ischemic stroke etiology and clot origin. \n* However, some difficulties arise/exist that create challenges you could lend a hand in solving:\n * unique data formats\n * image file sizes\n * the number of available pathology slides\n \n**S**troke **T**hromboembolism **R**egistry of **I**maging and **P**athology (**STRIP**) is a uniquely large multicenter project led by Mayo Clinic Neurovascular Lab with the aim of histopathologic characterization of thromboemboli of various etiologies and examining clot composition and its relation to mechanical thrombectomy revascularization.\n\n
COMPETITION HOST INFORMATION\n\nThe Mayo Clinic is a nonprofit American academic medical center focused on integrated health care, education, and research.\n\nTo decrease the chances of subsequent strokes, the Mayo Clinic Neurovascular Research Laboratory encourages data scientists to improve artificial intelligence-based etiology classification so that physicians are better equipped to prescribe the correct treatment. New computational and artificial intelligence approaches could help save the lives of stroke survivors and help us better understand the world's second-leading cause of death.\n\n
VISUAL EXPLANATION\n\n

\n\nFigure 1:   Image Showing an Acute Ischemic Stroke Clot
        **Left Image:** Example of an acute ischemic stroke clot stained with Martius scarlet blue stain;
        **Right Image:** Confocal image demonstrating the presence of platelets (red), fibrinogen (green) and white blood cells (blue).\n
\n\n\n
DETAILS OF Stroke Thromboembolism Registry of Imaging and Pathology\n\nWith the arrival of mechanical thrombectomy for the treatment of acute ischemic stroke, there has been growing interest in the histopathological and imaging characteristics of retrieved thromboemboli. Past studies propose that studying clot composition can provide insights into stroke etiology, as well as demonstrate association recanalization success with intravenous thrombolysis and mechanical thrombectomy. Furthermore, improved characterization of clot composition regarding imaging and histopathologic features could help in device selection prior to the mechanical thrombectomy procedure itself.\n\n**STRIP** is a multicenter project that aims to:\n* Characterize the histopathologic features of thromboemboli of various etiologies using a combination of histology (hematoxylin and eosin and Martius scarlet blue) and immunostaining methods\n* Examine the correlation between clot composition and degree and ease of revascularization using mechanical thrombectomy techniques\n* Review the association between device selection (aspiration versus stent retriever), technique selection (balloon guide catheter versus standard catheter), clot subtype (RBC rich, platelet rich, fibrin rich, WBC rich) and angiographic outcomes\n* Study the association between the presence of a hyperdense artery sign (on imaging) and clot histological characteristics\n\n
COMPETITION IMPACT STATEMENT\n\nStroke remains the second-leading cause of death worldwide. Each year in the United States, over 700,000 individuals experience an ischemic stroke caused by a blood clot blocking an artery to the brain. A second stroke (23% of total events are recurrent) worsens the chances of the patient’s survival. However, subsequent strokes may be mitigated if physicians can determine stroke etiology, which influences the therapeutic management following stroke events.\n\nYour work will enable healthcare providers to better identify the origins of blood clots in deadly strokes, making it easier for physicians to prescribe the best post-stroke therapeutic management and reducing the likelihood of a second stroke.","metadata":{}},{"cell_type":"markdown","source":"
\n\n

1.2 COMPETITION EVALUATION

\n\n---\n\n
GENERAL EVALUATION INFORMATION\n\nEach image has been labeled with an etiology class and for each image you must submit a **probability** for all classes.\n* **CE** – *CardioEmbolic*\n* **LAA** – *Large-Artery Atheroscle*\n\nSubmissions are evaluated using a weighted multi-class logarithmic loss. The overall effect is such that each class is roughly equally important for the final score.\n\n
\n\nThe formula is then:\n\n$$\n\\text{Log Loss} = - \\left( \\frac{\\sum^{M}_{i=1} w_{i} \\cdot \\sum_{j=1}^{N_{i}} \\frac{y_{ij}}{N_{i}} \\cdot \\ln p_{ij} }{\\sum^{M}_{i=1} w_{i}} \\right)\n$$\n\nWhere: \n* $N$ is the number of images in the class set\n* $M$ is the number of classes\n* $ln$ is the natural logarithm\n* $y_{ij}$ is $1$ if observation $i$ belongs to class $j$ and $0$ otherwise...\n* $p_{ij}$ is the predicted probability that image $i$ belongs to class $j$\n\n

\n\nThe submitted probabilities for a given image are not required to sum to one because they are rescaled prior to being scored (each row is divided by the row sum). In order to avoid the extremes of the log function, each predicted probability $p$ is replaced with $\\max(\\min(p,1-10^{-15}),10^{-15})$.\n\n
\n\n---\n\n
\n\n
SUBMISSION FILE INFORMATION\n\nFor each **`patient_id`** in the test set, you must predict a probability for each of the two etiology classes. \n\nThe file should contain a header and have the following format:\n\n```\npatient_id,CE,LAA\n01f2b3,0.5,0.5\n04de22,0.4,0.6\n0a47c9,0.99,0.01\n0af8b6,0.9999999,0.00000001\n...\n```\n\n\n
IS THIS A CODE COMPETITION?\n\nYES. AND YOU ONLY HAVE 1 SUBMISSION PER DAY!\n\n
","metadata":{}},{"cell_type":"markdown","source":"
\n\n

1.3 DATASET OVERVIEW

\n\n---\n\n
GENERAL INFORMATION\n\nThe dataset for this competition:\n* Over a thousand high-resolution whole-slide digital pathology images\n* Each slide depicts a blood clot from a patient that had experienced an acute ischemic stroke\n* The slides comprising the training and test sets depict clots with an etiology (that is, origin)\n * Either **CE** (Cardioembolic) \n * or **LAA** (Large Artery Atherosclerosis). \n* The hosts include a set of supplemental slides with either an unknown etiology or an etiology other than CE or LAA.\n\n
\n\n**Your task is to classify the etiology (CE or LAA) of the slides in the test set for each patient.**\n\n
FILE & DATA FIELD INFORMATION\n\n**`train/`** \n- A folder containing images in the TIFF format to be used as training data.\n\n**`test/`** \n- A folder containing images to be used as test data. \n- The actual test data (swapped in on submission) comprises about 280 images\n\n**`other/`** \n- A supplemental set of images with a either an unknown etiology or an etiology other than CE or LAA\n\n
\n\n**`train.csv`** \n- Contains annotations for images in the **`train/`** folder.\n- **Columns**\n - **`image_id`** \n - A unique identifier for this instance having the form {patient_id}_{image_num}. \n - Corresponds to the image `{image_id}.tif`.\n - **`center_id`**`** \n - Identifies the medical center where the slide was obtained.\n - **`patient_id `**\n - Identifies the patient from whom the slide was obtained.\n - **`image_num`** \n - Enumerates images of clots obtained from the same patient.\n - **`label`** \n - The etiology of the clot, either CE or LAA. \n - This field is the classification target.\n \n**`test.csv`**\n- Annotations for images in the **`test/`** folder. Has the same fields as **`train.csv`** excluding label.\n\n**`other.csv`**\n- Annotations for images in the **`other/`** folder. Has the same fields as **`train.csv`**. The `center_id` is unavailable for these images however.\n- **Columns**\n - **`label`** \n - The etiology of the clot, either `Unknown` or `Other`.\n - **`other_specified`** \n - The specific etiology, when known, in case the etiology is labeled as `Other`.\n\n**`sample_submission.csv`**\n- A sample submission file in the correct format. \n- Note in particular that you should make one prediction per **`patient_id`**, not per **`image_id`**.\n\n
","metadata":{}},{"cell_type":"markdown","source":"
\n\n\n\n

2  SETUP    ⤒

\n\n---\n","metadata":{}},{"cell_type":"markdown","source":"
\n\n

2.1 ACCELERATOR DETECTION

\n\n---\n\nIn order to use **`TPU`**, we use **`TPUClusterResolver`** for the initialization which is necessary to connect to the remote cluster and initialize cloud TPUs. Let's go over two important points\n\n1. When using TPU on Kaggle, you don't need to specify arguments for **`TPUClusterResolver`**\n2. However, on **G**oogle **C**ompute **E**ngine (**GCE**), you will need to do the following:\n\n
\n\n```python\n# The name you gave to the TPU to use\nTPU_WORKER = 'my-tpu-name'\n\n# or you can also specify the grpc path directly\n# TPU_WORKER = 'grpc://xxx.xxx.xxx.xxx:8470'\n\n# The zone you chose when you created the TPU to use on GCP.\nZONE = 'us-east1-b'\n\n# The name of the GCP project where you created the TPU to use on GCP.\nPROJECT = 'my-tpu-project'\n\ntpu = tf.distribute.cluster_resolver.TPUClusterResolver(tpu=TPU_WORKER, zone=ZONE, project=PROJECT)\n```\n\n
\n πŸ›‘   WARNING:

- Although the Tensorflow documentation says it is the project name that should be provided for the argument `project`, it is actually the Project ID, that you should provide. This can be found on the GCP project dashboard page.
\n
\n\n
\n πŸ“–   REFERENCES:

\n - Guide - Use TPUs
\n - Doc - TPUClusterResolver
\n\n
","metadata":{}},{"cell_type":"code","source":"print(f\"\\n... ACCELERATOR SETUP STARTING ...\\n\")\n\n# Detect hardware, return appropriate distribution strategy\ntry:\n # TPU detection. No parameters necessary if TPU_NAME environment variable is set. On Kaggle this is always the case.\n TPU = tf.distribute.cluster_resolver.TPUClusterResolver() \nexcept ValueError:\n TPU = None\n\nif TPU:\n print(f\"\\n... RUNNING ON TPU - {TPU.master()}...\")\n tf.config.experimental_connect_to_cluster(TPU)\n tf.tpu.experimental.initialize_tpu_system(TPU)\n strategy = tf.distribute.experimental.TPUStrategy(TPU)\nelse:\n print(f\"\\n... RUNNING ON CPU/GPU ...\")\n # Yield the default distribution strategy in Tensorflow\n # --> Works on CPU and single GPU.\n strategy = tf.distribute.get_strategy() \n\n# What Is a Replica?\n# --> A single Cloud TPU device consists of FOUR chips, each of which has TWO TPU cores. \n# --> Therefore, for efficient utilization of Cloud TPU, a program should make use of each of the EIGHT (4x2) cores. \n# --> Each replica is essentially a copy of the training graph that is run on each core and \n# trains a mini-batch containing 1/8th of the overall batch size\nN_REPLICAS = strategy.num_replicas_in_sync\n \nprint(f\"... # OF REPLICAS: {N_REPLICAS} ...\\n\")\n\nprint(f\"\\n... ACCELERATOR SETUP COMPLTED ...\\n\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n\n

2.2 COMPETITION DATA ACCESS

\n\n---\n\nTPUs read data must be read directly from **G**oogle **C**loud **S**torage **(GCS)**. Kaggle provides a utility library – **`KaggleDatasets`** – which has a utility function **`.get_gcs_path`** that will allow us to access the location of our input datasets within **GCS**.

\n\n
\n πŸ“Œ   TIPS:

- If you have multiple datasets attached to the notebook, you should pass the name of a specific dataset to the `get_gcs_path()` function. In our case, the name of the dataset is the name of the directory the dataset is mounted within.

\n
","metadata":{}},{"cell_type":"code","source":"print(\"\\n... DATA ACCESS SETUP STARTED ...\\n\")\n\nif TPU:\n # Google Cloud Dataset path to training and validation images\n DATA_DIR = KaggleDatasets().get_gcs_path('mayo-clinic-strip-ai')\n save_locally = tf.saved_model.SaveOptions(experimental_io_device='/job:localhost')\n load_locally = tf.saved_model.LoadOptions(experimental_io_device='/job:localhost')\nelse:\n # Local path to training and validation images\n DATA_DIR = \"/kaggle/input/mayo-clinic-strip-ai\"\n save_locally = None\n load_locally = None\n\nprint(f\"\\n... DATA DIRECTORY PATH IS:\\n\\t--> {DATA_DIR}\")\n\nprint(f\"\\n... IMMEDIATE CONTENTS OF DATA DIRECTORY IS:\")\nfor file in tf.io.gfile.glob(os.path.join(DATA_DIR, \"*\")): print(f\"\\t--> {file}\")\n\nprint(\"\\n\\n... DATA ACCESS SETUP COMPLETED ...\\n\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n\n

2.3 LEVERAGING XLA OPTIMIZATIONS

\n\n---\n\n\n**XLA** (Accelerated Linear Algebra) is a domain-specific compiler for linear algebra that can accelerate TensorFlow models with potentially no source code changes. **The results are improvements in speed and memory usage**.\n\n
\n\nWhen a TensorFlow program is run, all of the operations are executed individually by the TensorFlow executor. Each TensorFlow operation has a precompiled GPU/TPU kernel implementation that the executor dispatches to.\n\nXLA provides us with an alternative mode of running models: it compiles the TensorFlow graph into a sequence of computation kernels generated specifically for the given model. Because these kernels are unique to the model, they can exploit model-specific information for optimization.

\n\n
\n πŸ›‘   WARNING:

- XLA can not currently compile functions where dimensions are not inferrable: that is, if it's not possible to infer the dimensions of all tensors without running the entire computation
\n
\n\n
\n πŸ“Œ   NOTE:

- XLA compilation is only applied to code that is compiled into a graph (in TF2 that's only a code inside tf.function).
- The jit_compile API has must-compile semantics, i.e. either the entire function is compiled with XLA, or an errors.InvalidArgumentError exception is thrown)\n
\n\n
\n πŸ“–   REFERENCE:

- XLA: Optimizing Compiler for Machine Learning
\n
","metadata":{}},{"cell_type":"code","source":"print(f\"\\n... XLA OPTIMIZATIONS STARTING ...\\n\")\n\nprint(f\"\\n... CONFIGURE JIT (JUST IN TIME) COMPILATION ...\\n\")\n# enable XLA optmizations (10% speedup when using @tf.function calls)\ntf.config.optimizer.set_jit(True)\n\nprint(f\"\\n... XLA OPTIMIZATIONS COMPLETED ...\\n\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n\n

2.4 BASIC DATA DEFINITIONS & INITIALIZATIONS

\n\n---\n","metadata":{}},{"cell_type":"code","source":"print(\"\\n... BASIC DATA SETUP STARTING ...\\n\\n\")\n\n# Open the training dataframe and display the initial dataframe\nTRAIN_DIR = os.path.join(DATA_DIR, \"train\")\nTRAIN_CSV = os.path.join(DATA_DIR, \"train.csv\")\ntrain_df = pd.read_csv(TRAIN_CSV)\ntrain_df[\"image_path\"] = train_df[\"image_id\"].apply(lambda x: os.path.join(TRAIN_DIR, x+\".tif\"))\n\n# Add size information to train dataframe\ntrain_df[\"image_size\"] = train_df.image_path.progress_apply(lambda x: Image.open(x).size)\ntrain_df[\"image_width\"] = train_df[\"image_size\"].apply(lambda x: int(x[0]))\ntrain_df[\"image_height\"] = train_df[\"image_size\"].apply(lambda x: int(x[1]))\ntrain_df[\"image_area\"] = train_df[\"image_width\"]*train_df[\"image_height\"]\n\nprint(\"\\n... TRAINING DATAFRAME... \\n\")\ndisplay(train_df)\n\n# Open the testing dataframe and display the initial dataframe\nTEST_DIR = os.path.join(DATA_DIR, \"test\")\nTEST_CSV = os.path.join(DATA_DIR, \"test.csv\")\ntest_df = pd.read_csv(TEST_CSV)\ntest_df[\"image_path\"] = test_df[\"image_id\"].apply(lambda x: os.path.join(TEST_DIR, x+\".tif\"))\n\n# Add size information to test dataframe\ntest_df[\"image_size\"] = test_df.image_path.progress_apply(lambda x: Image.open(x).size)\ntest_df[\"image_width\"] = test_df[\"image_size\"].apply(lambda x: int(x[0]))\ntest_df[\"image_height\"] = test_df[\"image_size\"].apply(lambda x: int(x[1]))\ntest_df[\"image_area\"] = test_df[\"image_width\"]*test_df[\"image_height\"]\n\nprint(\"\\n... TESTING DATAFRAME... \\n\")\ndisplay(test_df)\n\n# Open the other dataframe and display the initial dataframe\nOTHER_DIR = os.path.join(DATA_DIR, \"other\")\nOTHER_CSV = os.path.join(DATA_DIR, \"other.csv\")\nother_df = pd.read_csv(OTHER_CSV)\n\nprint(\"\\n... OTHER DATAFRAME... \\n\")\ndisplay(other_df)\n\n# Open the sample submission dataframe\nSS_CSV = os.path.join(DATA_DIR, \"sample_submission.csv\")\nss_df = pd.read_csv(SS_CSV)\n\nprint(\"\\n... SAMPLE SUBMISSION DATAFRAME... \\n\")\ndisplay(ss_df)\n\n# For debugging purposes when the test set hasn't been substituted we will know\nDEBUG=len(ss_df)==4\nprint(f\"\\n\\n\\n... ARE WE DEBUGGING: {DEBUG}... \\n\")\n\n# Capture examples that span smallest, median, and large image size examples\nSMALLEST_IMAGE_ID = train_df[train_df.image_area==train_df.image_area.min()].image_id.values[0]\nSMALLEST_IMAGE_ROW = train_df[train_df.image_id==SMALLEST_IMAGE_ID]\nMEDIAN_IMAGE_ID = train_df.sort_values(by=\"image_area\").reset_index().iloc[len(train_df)//2].image_id\nMEDIAN_IMAGE_ROW = train_df[train_df.image_id==MEDIAN_IMAGE_ID]\nLARGEST_IMAGE_ID = train_df[train_df.image_area==train_df.image_area.max()].image_id.values[0]\nLARGEST_IMAGE_ROW = train_df[train_df.image_id==LARGEST_IMAGE_ID]\n\nprint(\"\\n... BASIC DATA SETUP FINISHED ...\\n\\n\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n\n\n\n\n\n

\n 3  HELPER FUNCTION & CLASSES    ⤒\n

\n\n---\n\nI've created some helper functions specifically for this competition to deal with the more unwieldy images. I tried my best to document everything thoroughly but if you need help please leave a comment!","metadata":{}},{"cell_type":"code","source":"# Different interpolcasefoldation based on version of Kaggle Kernel\ntry:\n DEFAULT_INTERPOLATION = Image.resampling.BILINEAR\nexcept Exception as AttributeError:\n DEFAULT_INTERPOLATION = Image.BILINEAR\n \n\ndef flatten_l_o_l(nested_list):\n \"\"\" Flatten a list of lists \"\"\"\n return [item for sublist in nested_list for item in sublist]\n\ndef pil_open_downscaled_slide(img_path, downsample_by=8, interpolation=None, reducing_gap=3.0, as_numpy=True):\n \"\"\"\n \n Helper function to convert WSI into smaller downscaled version using PIL\n \n Timing details for MAYO CLINIC STRIP AI dataset:\n SMALLEST IMAGE BY AREA (4417, 5314)\n * Function takes ~0.5 seconds to run\n MEDIAN IMAGE BY AREA (17573, 38743)\n * Function takes ~16 seconds to run\n LARGEST IMAGE BY AREA (48282, 101406)(~208X LARGER THAN SMALLEST IMAGE)(~7X LARGER THAN MEDIAN IMAGE)\n * Function will CRASH!\n \n Args:\n img_path (str): Path to .tif file to be downsampled\n downsample_by (int): How many times smaller should resultant \n image be. i.e. image_size*(1/downsample_by) = new_size\n interpolation (int, optional): The value that enumerates\n the type of interpolation to be used in downsampling.\n This can be one\n * `Image.NEAREST` --> ENUM=0\n * `Image.BOX` --> ENUM=1\n * `Image.BILINEAR` --> ENUM=2\n * `Image.HAMMING` --> ENUM=3\n * `Image.BICUBIC` --> ENUM=4\n * `Image.LANCZOS` --> ENUM=5\n If omitted, it defaults to `Image.BILINEAR`\n reducing_gap (float, optional): How many steps to use when interpolating/resampling\n the image after integer downscaling is performed. \n ELI5: \n * Higher (Approaching 5) means it takes longer but is better quality\n * Lower (Approaching 0) means it is much faster but lower quality\n as_numpy (bool, optional): Whether to return image as numpy array (default)\n or leave as PIL.Image object for further manipulation\n \n Returns:\n Downsampled image as a numpy array of type uint8 with only 3 channels\n \"\"\"\n \n # Catches and warning\n if downsample_by<8: print(\"\\n... WARNING – DUE TO LOW DOWNSCALE_BY VALUE THE RETURND IMAGE MAY OCCUPY A LARGE AMOUNT OF RAM – WARNING ...\\n\")\n if interpolation is None: interpolation=DEFAULT_INTERPOLATION\n \n # Open the image with PIL\n tmp_img = Image.open(img_path)\n \n # Get original dimensions and calculate new dimensions\n original_dimensions = tmp_img.size\n downsample_dimensions = (int(original_dimensions[0]/downsample_by), int(original_dimensions[1]/downsample_by))\n \n # Resize the image\n tmp_img.thumbnail(size=downsample_dimensions, \n resample=interpolation,\n reducing_gap=reducing_gap)\n \n # Return either np.ndarray or PIL.Image object\n return np.asarray(tmp_img) if as_numpy else tmp_img\n \ndef pyvips_open_downsampled_slide(img_path, downsample_by=8, as_numpy=True, verbose=True):\n \"\"\"\n \n Helper function to convert WSI into smaller downscaled version using pyvips.\n \n Timing details for MAYO CLINIC STRIP AI dataset:\n SMALLEST IMAGE BY AREA (4417, 5314)\n * Function takes ~1 seconds to run\n MEDIAN IMAGE BY AREA (17573, 38743)(~30X LARGER THAN SMALLEST IMAGE)\n * Function takes ~30 seconds to run\n LARGEST IMAGE BY AREA (48282, 101406)(~208X LARGER THAN SMALLEST IMAGE)(~7X LARGER THAN MEDIAN IMAGE)\n * Function takes ~405 seconds to run\n \n Args:\n img_path (str): Path to .tif file to be downsampled\n downsample_by (int): How many times smaller should resultant \n image be. i.e. image_size*(1/downsample_by) = new_size\n as_numpy (bool, optional): Whether to return image as numpy array (default)\n or leave as PIL.Image object for further manipulation\n \n Returns:\n Downsampled image as a numpy array of type uint8 with only 3 channels\n \"\"\"\n \n # Open the image with PIL\n tmp_img = pyvips.Image.new_from_file(img_path) \n \n if verbose:\n print(\"\\n... APPROXIMATE TIME TO LOAD IMAGE IS AT MOST APPROXIMATELY: \" \\\n f\"{int((405/(48282*101406))*(tmp_img.width*tmp_img.height))} SECONDS ...\\n\")\n \n # Resize the image\n tmp_img = tmp_img.resize(1/downsample_by)\n \n # Return either np.ndarray or PIL.Image object\n return tmp_img.numpy() if as_numpy else tmp_img\n \ndef get_slide_tiles(img_path, tile_size=(384,384), stride=(256, 256), downsample_by=4, drop_empty=True, drop_near_empty=True, near_empty_sum_thresh=2500, near_empty_bg_thresh=170, return_downsampled_image=False, discard_non_square=True, verbose=True):\n \"\"\" Tile a WSI\n \n Args:\n img_path (str): Path to .tif file to be downsampled \n tile_size (tuple of ints, optional): The dimensions of the tile/patches to be extracted\n stride (tuple of ints, optional): The stride between tile/patches. \n i.e. a stride of two: xx____, __xx__, ____xx\n downsample_by (int, optional): How many times smaller should resultant \n image be. i.e. image_size*(1/downsample_by) = new_size\n drop_empty (bool, optional): Whether or not to drop images that are determined to be empty\n i.e. ALL the pixels in the tile are determined to be background pixels\n drop_near_empty (bool, optional): Whether or not to drop images that are determined to be \"mostly\" empty\n i.e. MOST of the pixels in the tile are determined to be background pixels\n near_empty_sum_thresh (int, optional): The threshold to use to determine what \"MOST\"\n means in the context of background vs forgeround.\n i.e. How many pixels are allowed to be background pixels before we drop?\n near_empty_bg_thresh (int, optional): What qualifies as the threshold at which if the pixel\n is brighter than this it is considered background... wheras lower than this is foreground\n (I AM AWARE FOREGROUND WILL YIELD SOME PERCENTAGE OF \"BACKGROUND\" PIXELS)\n return_downsampled_image (bool, optional): Whether we return the downsampled numpy array\n discard_non_square (bool, optional): Whether to discard or pad non-square edge tiles\n \n Returns:\n A list of tiles – np.ndarray of shape `tile_size`\n \"\"\"\n \n o_img = pyvips_open_downsampled_slide(img_path, downsample_by, verbose=verbose)\n (orig_h, orig_w), (tile_h, tile_w), (stride_w, stride_h) = o_img.shape[:-1], tile_size, stride\n \n tile_map, non_bg_sum = {}, near_empty_sum_thresh\n \n row_steps = range(0, orig_h+stride_h, stride_h)\n col_steps = range(0, orig_w+stride_w, stride_w)\n n_steps = len(row_steps)*len(col_steps)\n for j, y in enumerate(row_steps):\n for i, x in enumerate(col_steps):\n tile = o_img[y:(y+tile_h), x:(x+tile_w)]\n if drop_empty or drop_near_empty: \n non_bg_sum = np.where(tile>near_empty_bg_thresh, 0, 1).sum()\n if (drop_empty and non_bg_sum==0) or (drop_near_empty and non_bg_sum03} \"\n if not discard_non_square or tile.shape[0]==tile.shape[1]:\n tile_map[loc_str]=tile\n if (n_steps<=999 and verbose): print(loc_str, end=\"\")\n else:\n if (n_steps<=999 and verbose): print(\" XXX \", end=\"\")\n if (n_steps<=999 and verbose): print(\"\\n\")\n \n if return_downsampled_image:\n return tile_map, o_img\n else:\n return tile_map\n\ndef process_and_plot(pd_row, draw_grid=True, grid_stride=(256,256), grid_size=(384,384), draw_text=True):\n def _add_text(img, annot, color, pt1, pt2, **kwargs):\n text_width, text_height = cv2.getTextSize(annot, FONT, FONT_SCALE, FONT_THICKNESS)[0]\n annot_loc = (int(round(pt1[0]+((pt2[0]-pt1[0])/2)-(text_width/2))), int(round(pt1[1]+((pt2[1]-pt1[1])/2)+(text_height/2))))\n img = cv2.putText(img, annot, annot_loc, FONT, FONT_SCALE, color, FONT_THICKNESS, FONT_LINE_TYPE)\n return img\n \n print(f\"\\n... PROCESSING IMAGE {pd_row.image_id.values[0]} - {pd_row.image_size.values[0]} ...\\n\")\n t1 = time.time()\n image_tile_map, downsampled_image = get_slide_tiles(pd_row.image_path.values[0], return_downsampled_image=True)\n \n n_tiles = len(image_tile_map)\n print(f\"\\nIMAGE PROCESSING ELAPSED TIME: {time.time()-t1:.4f} SECONDS\")\n print(\"RESULTANT NUMBER OF TILES: \", n_tiles)\n \n print(\"\\n\\n... ORIGINAL IMAGE - DOWNSAMPLED BY A FACTOR OF 4 ...\\n\\n\")\n w_to_h_ratio = pd_row.image_width.values[0]/pd_row.image_height.values[0]\n plt.figure(figsize=(int(20*w_to_h_ratio), 20) if w_to_h_ratio>1 else (20, int(20/w_to_h_ratio)))\n plt.axis(False)\n plt.title(f\"Complete Downsampled Image – {pd_row.image_id.values[0]}\", fontweight=\"bold\")\n if draw_grid:\n clr_cnt, c_map, clrs, draw_map_fg, draw_map_bg = 0, plt.get_cmap(\"rainbow\"), [], [], []\n row_steps = range(0, downsampled_image.shape[0]+grid_stride[1], grid_stride[1])\n col_steps = range(0, downsampled_image.shape[1]+grid_stride[0], grid_stride[0])\n clrs = [[int(x*255) for x in c_map((i+1)/len(image_tile_map.keys()))] for i in range(len(image_tile_map.keys()))]\n draw_image = downsampled_image.copy()\n if draw_text: \n draw_annots, FONT, FONT_SCALE, FONT_THICKNESS, FONT_LINE_TYPE = [], cv2.FONT_HERSHEY_SIMPLEX, 1.56789, 3, cv2.LINE_AA\n \n for j, y in enumerate(row_steps):\n for i, x in enumerate(col_steps):\n if f\" {j*len(col_steps)+i+1:>03} \" in image_tile_map.keys():\n draw_map_fg.append(dict(pt1=(x,y), pt2=(x+grid_size[0], y+grid_size[1]), color=clrs[clr_cnt], thickness=3, annot=f\" {j*len(col_steps)+i+1:>03} \"))\n clr_cnt +=1\n else:\n draw_map_bg.append(dict(pt1=(x,y), pt2=(x+grid_size[0], y+grid_size[1]), color=(10,10,10), thickness=2, annot=f\" {j*len(col_steps)+i+1:>03} \")) \n \n for draw_map in draw_map_bg+draw_map_fg:\n draw_image = cv2.rectangle(draw_image, draw_map[\"pt1\"], draw_map[\"pt2\"], draw_map[\"color\"], draw_map[\"thickness\"])\n if draw_text:\n _add_text(draw_image, draw_map[\"annot\"], draw_map[\"color\"], draw_map[\"pt1\"], draw_map[\"pt2\"])\n ### STRIDE LINES ###\n # hline_list = list(range(0, downsampled_image.shape[0]+sample_grid_stride[0], sample_grid_stride[0]))\n # vline_list = list(range(0, downsampled_image.shape[1]+sample_grid_stride[1], sample_grid_stride[1]))\n # plt.hlines(hline_list, xmin=0, xmax=sample_grid_stride[0]*(len(vline_list)-1), colors=\"k\")\n # plt.vlines(vline_list, ymin=0, ymax=sample_grid_stride[1]*(len(vline_list)-1), colors=\"k\")\n \n plt.imshow(downsampled_image if not draw_grid else draw_image)\n plt.tight_layout()\n plt.show()\n\n print(f\"\\n\\n\\n... PLOTTING FIRST {min(16, n_tiles)} OF {n_tiles} TILES FROM IMAGE ...\\n\")\n plt.figure(figsize=(20,20))\n for i, (loc_str, tile) in enumerate(image_tile_map.items()):\n if i==16: break\n plt.subplot(4,4,i+1)\n plt.title(f\"Tile Location: {loc_str}\", fontweight=\"bold\")\n plt.axis(False)\n plt.imshow(tile)\n plt.tight_layout()\n plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"\\n... TEST WITH SMALLEST IMAGE BY AREA ...\\n\")\nprocess_and_plot(SMALLEST_IMAGE_ROW)\n\nprint(\"\\n... TEST WITH MEDIAN IMAGE BY AREA ...\\n\")\nprocess_and_plot(MEDIAN_IMAGE_ROW)\n\nprint(\"\\n... TEST WITH LARGEST IMAGE BY AREA ...\\n\")\nprocess_and_plot(LARGEST_IMAGE_ROW)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n\n\n\n\n\n

\n 4  DATASET EXPLORATION    ⤒\n

\n\n---","metadata":{}},{"cell_type":"markdown","source":"
\n\n

4.0 LOOK AT A FEW EXAMPLES FOR EACH CLASS

\n\n---\n\nWe simply do this to make sure everything is where it should be and we understand the basics of how to access all the relevant data.","metadata":{}},{"cell_type":"code","source":"n_demo_examples = 5\ndemo_df = train_df.drop_duplicates(\"patient_id\").head(6).sort_values(by=[\"label\", \"image_area\"]).reset_index(drop=True)\nce_demo_rows = train_df[train_df.label==\"CE\"].drop_duplicates(\"patient_id\").sort_values(by=\"image_area\").head(n_demo_examples).reset_index(drop=True)\nlaa_demo_rows = train_df[train_df.label==\"LAA\"].drop_duplicates(\"patient_id\").sort_values(by=\"image_area\").head(n_demo_examples).reset_index(drop=True)\n\nprint(\"\\n\\n... CE DEMO ROWS ...\\n\\n\\n\")\ndisplay(ce_demo_rows)\n\nprint(\"\\n\\n... CE DEMO ROWS WHOLE IMAGE VISUALIZED AND TILES VISUALIZED ...\\n\")\nfor i in range(len(ce_demo_rows)):\n process_and_plot(ce_demo_rows.iloc[i:i+1])\n\nprint(\"\\n\\n\\n\\n... LAA DEMO ROWS ...\\n\\n\\n\")\ndisplay(laa_demo_rows)\n\nprint(\"\\n\\n... LAA DEMO ROWS WHOLE IMAGE VISUALIZED AND TILES VISUALIZED ...\\n\")\nfor i in range(len(laa_demo_rows)):\n process_and_plot(laa_demo_rows.iloc[i:i+1])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n\n

4.1 INVESTIGATE THE `IMAGE_ID` COLUMN

\n\n---\n\n**The `image_id` value is a unique identifier for this instance having the form `{patient_id}_{image_num}`.** This is simply a UUID associated with every example. Every row in every dataset train/test/public/private should have a unique value in for this field. \n\nIt can be very useful in identifying examples if we end up tiling or processing the dataset.","metadata":{}},{"cell_type":"code","source":"print(f\"\\n... THERE ARE {train_df.image_id.nunique()} UNIQUE VALUES FOR {len(train_df)} EXAMPLES IN THE TRAINING DATA ...\")\nprint(f\"... THERE ARE {test_df.image_id.nunique()} UNIQUE VALUES FOR {len(test_df)} EXAMPLES IN THE TEST DATA ...\\n\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n\n

4.2 INVESTIGATE THE `CENTER_ID` COLUMN

\n\n---\n\n**The `center_id` value identifies the medical center where the slide was obtained.** There are 11 different **`center_id`** values in the dataset. \n\n
OBSERVATIONS\n* A large number of slides was obtained from a single medical center:\n * **`center_id==257`**\n * 257 of the 754 examples were obtained here (34.1%)\n* No center has fewer than 16 examples associated with it (2.1%)","metadata":{}},{"cell_type":"code","source":"print(f\"\\n... `center_id` HAS {train_df['center_id'].nunique()} UNIQUE VALUES FOR {len(train_df)} EXAMPLES IN THE TRAINING DATA ...\\n\\n\")\ntrain_df[\"center_id_count\"] = train_df.groupby(\"center_id\")[\"image_id\"].transform('count')\nfig = px.histogram(train_df, \"center_id\", \n color=\"center_id_count\", \n color_discrete_sequence=px.colors.sequential.Inferno_r,\n title=\"Histogram Showing Frequency of `center_id` Values In Training Dataset\",\n labels={\"color\":\"Center ID Legend\"},\n )\nfig.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n\n

4.3 INVESTIGATE THE `PATIENT_ID` COLUMN

\n\n---\n\n**The `patient_id` value identifies the patient from whom the slide was obtained.** There are 632 different **`patient_id`** values in the dataset (that's almost 84% of the training examples!). \n\n
OBSERVATIONS\n* 543 of the patients are only represented by a single example (72.1%).\n* Globally, we can see that the **`center_id`** distributions are normally distributed across patients with 1 through 5 observations.","metadata":{}},{"cell_type":"code","source":"print(f\"\\n... `patient_id` HAS {train_df['patient_id'].nunique()} UNIQUE VALUES FOR {len(train_df)} EXAMPLES IN THE TRAINING DATA ...\\n\\n\")\ntrain_df[\"patient_id_count\"] = train_df.groupby(\"patient_id\")[\"image_id\"].transform('count')\nfig = px.histogram(train_df, \"patient_id_count\", color=\"center_id\", \n title=\"Histogram Showing Frequency of `patient_id` Values In Training Dataset Coloured By Center ID\",\n labels={\"center_id\":\"Center ID Legend\"})\nfig.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n\n

4.4 INVESTIGATE `IMAGE_NUM`

\n\n---\n\n**The `image_num` value enumerates images of clots obtained from the same patient.** This column is essentially repeated information along with **`patient_id`**. Unless there is some value in the enumeration, it simply allows us to build a unique identifier for a given patient's example.","metadata":{}},{"cell_type":"markdown","source":"
\n\n

4.5 INVESTIGATE `LABEL`

\n\n---\n\nThis is the big one!! **The `label` value gives the etiology of the clot, either CE or LAA. This field is the classification target.**\n\n
OBSERVATIONS\n\n* There are 2.64 times more annotations for CE than there are for LAA. Or in other words... 72.5% of the data is labelled as CE while only 27.5% of the data is labelled as LAA.\n* There doesn't appear to be a strong correlation between image size and label","metadata":{}},{"cell_type":"code","source":"fig = px.pie(train_df, \"label\", title=\"Pie Chart Showing Class Imbalance In Provided Label Clot Etiology\", hole=0.30901699, )\nfig.show()\n\nfig = px.histogram(train_df, \"image_area\", color=\"label\",\n title=\"Scatter Plot Showing Class Distribution According to Slide Size\")\nfig.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n\n

4.6 INVESTIGATE `IMAGE_SIZE`

\n\n---\n\nThis is a column created during the EDA, but the data is still relevant and important.\n\n
OBSERVATIONS\n\n* The image sizes are incredibly varied\n* The image sizes often trend square but on average are more vertical than horizontal\n* Image sizes range from $(4417, 5314)$ to $(48282, 101406)$... this is a MASSIVE difference with the larger image being nearly 209x times larger than the smaller image!!!!\n * To put it in perspective if we compared a 240p image vs a 2160p (4k) image... that's still less than a 100x increase... wheras in this dataset we have images that are 209x!","metadata":{}},{"cell_type":"code","source":"display(train_df[(train_df.image_area==train_df.image_area.min()) | (train_df.image_area==train_df.image_area.max())])\n\nfig = px.scatter(train_df, \"image_width\", \"image_height\", color=\"label\", size=\"image_area\",\n title=\"Scatter Plot Showing Class Distribution According to Slide Size\",\n width=1200, height=1200)\nfig.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**DATASET CREATION**","metadata":{}},{"cell_type":"code","source":"tiling_kwargs = dict(\n tile_size = (512,512),\n stride = (384,384),\n downsample_by = 6,\n drop_empty=True,\n drop_near_empty=True,\n near_empty_sum_thresh=50_000,\n near_empty_bg_thresh=170,\n return_downsampled_image=False,\n verbose=False\n)\n\ntile_image_ids, tile_patient_ids =[],[]\ntile_labels, tile_strs, tile_img_names =[],[],[]\nfor i, (_, row) in tqdm(enumerate(train_df.iterrows()), total=len(train_df)):\n tile_map = get_slide_tiles(row.image_path, **tiling_kwargs) \n for tile_str, tile_arr in tqdm(tile_map.items(), total=len(tile_map)):\n save_name = f\"{row['image_id']}_{row['patient_id']}__{tile_str.strip()}.png\"\n \n # Do saving (RGB to BGR because silly CV2)\n cv2.imwrite(save_name, tile_arr[..., ::-1])\n \n # Update dataframe\n tile_image_ids.append(row[\"image_id\"])\n tile_patient_ids.append(row[\"patient_id\"])\n tile_labels.append(row[\"label\"])\n tile_strs.append(tile_str.strip())\n tile_img_names.append(save_name)\ntile_df = pd.DataFrame({\"image_id\":tile_image_ids, \"patient_id\":tile_patient_ids, \"label\":tile_labels, \"tile_str\":tile_strs, \"img_name\":tile_img_names})","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}