{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":71549,"databundleVersionId":8561470,"sourceType":"competition"},{"sourceId":7457180,"sourceType":"datasetVersion","datasetId":4313900},{"sourceId":8650305,"sourceType":"datasetVersion","datasetId":5181481}],"dockerImageVersionId":30732,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# Run this to enable CSS types\nfrom IPython.core.display import HTML\n\n# Font stuff\nfont_to_use = \"Arsenal\" # \"Verdana\"\nfallback = \"Verdana\"\nfont_import_str = f\"\"\"\n@import url('https://fonts.googleapis.com/css2?family={font_to_use.replace(' ', '+')}:ital,wght@0,100;0,300;0,400;0,700;0,900;1,100;1,300;1,400;1,700;1,900&display=swap');\n\"\"\"\n\ndef css_styling(verbose=True):\n styles = open(\"/kaggle/input/my-css-styles/kaggle_styles.css\", \"r\").read().replace('Verdana', font_to_use) #+f\", {fallback}\")\n html_str = \"\"\n if verbose: print(html_str)\n return HTML(html_str)\n\ncss_styling(False)","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-06-23T19:01:09.447668Z","iopub.execute_input":"2024-06-23T19:01:09.448381Z","iopub.status.idle":"2024-06-23T19:01:09.470478Z","shell.execute_reply.started":"2024-06-23T19:01:09.448346Z","shell.execute_reply":"2024-06-23T19:01:09.469202Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n\n
\n\n
\n\n

CREATED BY: DARIEN SCHETTLER

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

THIS IS A WORK IN PROGRESS
\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\n
⚠️ WIP - UNDERGOING FREQUEST UPDATES - WIP ⚠️
\n\n
","metadata":{}},{"cell_type":"markdown","source":"

\n CHANGELOG\n

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

\n\n

\n TABLE OF CONTENTS\n

\n\n
\n\n

1    INTRODUCTION & JUSTIFICATION

\n\n
\n\n

2    BACKGROUND INFORMATION

\n\n
\n\n

3    IMPORTS

\n\n
\n\n

4    SETUP AND HELPER FUNCTIONS

\n\n
\n\n

5    EXPLORATORY DATA ANALYSIS

\n\n
\n\n

6    BASELINE SUBMISSION

\n\n
\n\n

7    CROSS VALIDATION

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

1  INTRODUCTION & JUSTIFICATION    ⤒

","metadata":{}},{"cell_type":"markdown","source":"

1.1 WHAT IS THIS?

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

1.2 WHY IS THIS?

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

1.3 WHO IS THIS FOR?

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

1.4 HOW WILL THIS WORK?

\n
\n\n\n

I'm going to assemble some markdown cells (like this one) at the beginning of the notebook to go over some concepts/details/etc.

\n\n

Following this, I will attempt to walk through the data and understand it better prior to composing a baseline solution.

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

2  BACKGROUND INFORMATION    ⤒

\n\n
\n\nLow back pain is a significant global health issue, impacting 619 million individuals in 2020, as reported by the World Health Organization. It's a common ailment that most people experience at some stage in their lives, with its prevalence increasing with age. Spondylosis, a group of degenerative spine conditions, often manifests through pain and restricted mobility. \n\nThe Spondylosis condition involves:\n- Degeneration of intervertebral discs\n- Narrowing of the spinal canal (spinal stenosis)\n- Narrowing of subarticular recesses or neural foramina\n- Compression or irritation of nerves in the lower back\n\nMagnetic resonance imaging (MRI) plays a crucial role in diagnosing these conditions, offering detailed views of the lumbar spine's vertebrae, discs, and nerves. Accurate diagnosis and grading of these conditions are essential for:\n- Guiding treatment options\n- Potential surgical interventions\n- Alleviating back pain\n- Enhancing patients' overall health and quality of life\n\nIn an effort to explore the potential of artificial intelligence in this field, the **Radiological Society of North America (RSNA)** has partnered with the **American Society of Neuroradiology (ASNR)** to launch **this competition**. \n\n
","metadata":{}},{"cell_type":"markdown","source":"

2.1 COMPETITION OVERVIEW

\n
\n\n
\n\nPRIMARY TASK DESCRIPTION\n
\n
\nThe goal of this competition is to create models that can be used to aid in the detection and classification of degenerative spine conditions using lumbar spine MR images. \n\nCompetitors will develop models that simulate a radiologist's performance in diagnosing spine conditions.\n\n
\n\nHOST TASK DESCRIPTION\n
\n
\nThis challenge aims to determine whether we can effectively detect and classify degenerative spine conditions using MRI images of the lumbar spine.\n\nThe focus of the challenge is on classifiying five specific lumbar spine degenerative conditions:\n1. Left Neural Foraminal Narrowing\n2. Right Neural Foraminal Narrowing\n3. Left Subarticular Stenosis\n4. Right Subarticular Stenosis\n5. Spinal Canal Stenosis\n\nFor each imaging study in the dataset, severity scores (Normal/Mild, Moderate, or Severe) are provided for each of these five conditions across the intervertebral disc levels:\n- L1/L2\n- L2/L3\n- L3/L4\n- L4/L5\n- L5/S1\n\nTo establish a robust ground truth dataset, the **RSNA challenge planning task force** collected imaging data from **eight sites across five continents**! This multi-institutional, expertly curated dataset aims to:\n- Standardize the classification of degenerative lumbar spine conditions\n- Facilitate the development of tools for rapid and accurate disease classification\n\n
\n\nVISUAL TASK DESCRIPTION AND ADDITIONAL BACKGROUND - SOURCED FROM THIS NOTEBOOK W/ ADDITIONAL DETAILS\n
\n
\nANATOMICAL OVERVIEW\n\nThe spine is a complex structure divided into four regions:\n- **Cervical region**: Consists of 7 vertebral bodies, located in the neck.\n- **Thoracic region**: Contains 12 vertebral bodies, located in the upper and mid-back.\n- **Lumbar region**: Comprises 5 vertebral bodies, found in the lower back.\n- **Sacral region**: Made up of 3-5 fused vertebral bodies, situated at the base of the spine.\n\n\"Spine\nFaculty of Washington University

\n\nEach region, except the sacrum, has **intervertebral discs** between the vertebral bodies. \n* These discs are cartilage structures that serve as cushions, absorbing shock and facilitating movement. \n\nRunning along the posterior aspect (back side) of each vertebral body is the **spinal cord**.\n* The **spinal cord** is a bundle of nerves that transmits signals between the brain and the rest of the body. \n* **Spinal nerves** exit the spinal cord through openings called **foramina**.\n\n\"Herniated\"\nMiami Neuroscience Center

\n\nCompression of the spinal cord or spinal nerves can cause significant pain. Factors contributing to this compression include:\n* **Bulging or herniated discs**: \n * When the intervertebral discs protrude out of their normal space.\n* **Degenerative changes in vertebrae**: \n * Leading to bony protrusions or bone spurs (osteophytes).\n* **Trauma**: \n * Physical injury that impacts the spine.\n* **Thickening of surrounding ligaments**: \n * The bands of tissue that connect bones and support the spine can thicken and compress nerves.\n\n
\n\nFORMINAL NARROWING OVERVIEW\n\nForamina, the openings through which spinal nerves exit, can become narrowed, a condition best visualized in the sagittal plane (a vertical slice parallel to the spine). \n* Foraminal narrowing results in nerve compression, leading to pain along the affected nerve's distribution, which is the path the nerve travels in the body.\n\n
\n
\n \"Foraminal\n
\n
\n \"Grading\n
\n


\n\nSUBARTICULAR STENOSIS OVERVIEW\n\n**Subarticular stenosis** occurs when there is compression of the spinal cord in the subarticular zone, a region just below the articular surface of the vertebrae. This condition is best visualized in the axial plane (a horizontal slice perpendicular to the spine). The compression in this zone is often due to similar factors that cause foraminal narrowing.\n\n
\n
\n \"Subarticular\n
\n
\n \"Grading\n
\n
\n\nLeft image from Miami Neuroscience Center

\n\nCANAL STENOSIS OVERVIEW\n\n**Canal stenosis** involves impingement of the spinal canal, the passageway through which the spinal cord travels. Causes include bulging discs, trauma, bony outgrowths (osteophytes), and thickened ligaments. The degree of compression is typically assessed in the axial plane.\n\n
\n
\n \"Canal\n
\n
\n \"Grading\n
\n
\n\nFrom Radiopedia

\n\nIMAGING OVERVIEW\n\n**MRI imaging** (Magnetic Resonance Imaging) of the spine can be performed in three planes: axial, sagittal, and coronal. For this challenge, the axial and sagittal planes are most relevant:\n\n- **Axial plane**: Horizontal slices perpendicular to the spine.\n- **Sagittal plane**: Vertical slices parallel to the spine.\n\nMRI images come in two primary types:\n\n- **T1-weighted images**: \n - These images highlight fat, making the inner parts of bones appear brighter.\n- **T2-weighted images**: \n - These images highlight water, making the spinal canal appear brighter.\n\n\nCase Western Reserve University

\n\nUnlike CT images (Computed Tomography), MRI images lack standardized pixel values, requiring careful consideration for standardization during analysis.\n\n

","metadata":{}},{"cell_type":"markdown","source":"

2.2 DATASET OVERVIEW

\n
\n\n
\n\n
HIGH LEVEL DATA SUMMARY\n\nThis is a complex challenge, requiring participants to predict the presence and severity of a set of degenerative conditions at each vertebral level of the lumbar spine. \n\nThe competition dataset includes more than 2,500 magnetic resonance imaging studies from eight sites on five continents. It has been annotated by a group of more than 50 expert radiologists recruited by the American Society of Neuroradiology.\n\nTo create the ground truth dataset, the RSNA challenge planning task force collected imaging data sourced from eight sites on five continents. This multi-institutional, expertly curated dataset promises to improve standardized classification of degenerative lumbar spine conditions and enable development of tools to automate accurate and rapid disease classification.\n\nThis competition uses a hidden test. When your submitted notebook is scored, the actual test data (including a full length sample submission) will be made available to your notebook.\n\n
\n\n
DATA COMPOSITION\n\nThe training dataset consists of 1975 rows (each row is a study)\n* Each study may include multiple series of images\n* Each study maps to target labels, such as `spinal_canal_stenosis_l1_l2`, with the severity levels of Normal/Mild, Moderate, or Severe. Some entries have incomplete labels.\n\nThe test dataset contains approximately 500 rows/studies (speculative based on host saying the dataset was ~2500 studies)\n\n
\n\n
DIRECTORY STRUCTURE\n\n```\nβ”œβ”€β”€ /kaggle/working\n└── /kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/\n β”œβ”€β”€ test_images/\n β”‚ β”œβ”€β”€ 1005139/\n β”‚ β”‚ └── 609308237/\n β”‚ β”‚ β”œβ”€β”€ 1.dcm\n β”‚ β”‚ └── ...\n β”‚ └── ...\n β”œβ”€β”€ test_series_descriptions.csv\n β”œβ”€β”€ train_images/\n β”‚ β”œβ”€β”€ 4003253/\n β”‚ β”‚ └── 702807833/\n β”‚ β”‚ β”œβ”€β”€ 1.dcm\n β”‚ β”‚ └── ...\n β”‚ └── ...\n β”œβ”€β”€ train_label_coordinates.csv\n β”œβ”€β”€ train_series_descriptions.csv\n └── train.csv\n```\n\n
\n\n
DATA FILE DESCRIPTIONS\n\ntrain.csv:\n* Contains labels for the training set. Each record includes:\n - study_id (string): The study ID. \n - Each study may include multiple series of images.\n - [condition]_[level] (string): The target labels, such as... \n - **`spinal_canal_stenosis_l1_l2`**, with severity levels of **`Normal/Mild`**, **`Moderate`**, or **`Severe`**. \n - Some entries have incomplete labels.\n\ntrain_label_coordinates.csv:\n* Provides the coordinates for labeled regions. Each record includes:\n - study_id (string): The study ID.\n - series_id (string): The imagery series ID.\n - instance_number (int): The image's order number within the 3D stack.\n - condition (string): The core condition, which can be one of...\n - **`spinal_canal_stenosis`**\n - **`neural_foraminal_narrowing`** [considered for each side of the spine]\n - **`subarticular_stenosis`** [considered for each side of the spine]\n - level (string): The relevant vertebrae, such as **`l3_l4`**.\n - x (float): The x-coordinate for the center of the labeled area.\n - y (float): The y-coordinate for the center of the labeled area.\n\nsample_submission.csv:\n* Provides a format template for submissions. Each record includes:\n - row_id (string): A slug of the study ID, condition, and level, such as `12345_spinal_canal_stenosis_l3_l4`.\n - normal_mild (float): The predicted probability for the Normal/Mild severity level.\n - moderate (float): The predicted probability for the Moderate severity level.\n - severe (float): The predicted probability for the Severe severity level.\n\n[train/test]_images/[study_id]/[series_id]/[instance_number].dcm:\n* The directory structure for the imagery data.\n\n[train/test]_series_descriptions.csv:\n* Contains descriptions of the scan series. Each record includes:\n - study_id (string): The study ID.\n - series_id (string): The series ID.\n - series_description (string): The scan's orientation.\n\n
\n\n
SAMPLE SUBMISSION EXAMPLE\n\n
row_id,normal_mild,moderate,severe\n123456_left_neural_foraminal_narrowing_l1_l2,0.333,0.333,0.333\n123456_left_neural_foraminal_narrowing_l2_l3,0.333,0.333,0.333\n123456_left_neural_foraminal_narrowing_l3_l4,0.333,0.333,0.333\netc.
\n\n
","metadata":{}},{"cell_type":"markdown","source":"

2.3 EVALUATION

\n
\n\n
\n\n
EVALUATION OVERVIEW\n\nSubmissions are evaluated using the average of sample weighted log losses and an any_severe_spinal prediction generated by the metric. The metric notebook can be found here.\n\nThe sample weights are as follows:\n- 1 for normal/mild.\n- 2 for moderate.\n- 4 for severe.\n\nFor each row ID in the test set, you must predict a probability for each of the different severity levels. \n\nNOTE: In rare cases the lowest vertebrae aren't visible in the imagery. You still need to make predictions (nulls will cause errors), but those rows will not be scored.\n\nFor this competition, the any_severe_scalar has been set to 1.0.\n\nThe submission file should contain a header and have the following format:\n\n
\nrow_id,normal_mild,moderate,severe\n123456_left_neural_foraminal_narrowing_l1_l2,0.333,0.333,0.333\n123456_left_neural_foraminal_narrowing_l2_l3,0.333,0.333,0.333\n123456_left_neural_foraminal_narrowing_l3_l4,0.333,0.333,0.333\netc.\n
\n\n
\n\n
\n\n
EVALUATION UPDATES (POST COMPETITION LAUNCH)\n\nMETRIC PATCH 1\n\nPost from Sohier Dane - Thread Here\n\n

As pointed out in this thread, the metric contained a bug that essentially eliminated the impact of the any_severe_spinal calculation. It was the same for all submissions. I've posted a preliminary patch to the metric and will rescore all submissions after I have time to conduct another review of the metric to make sure that was the only oversight.

\n\n
\n\nMETRIC PATCH 2\n\nPost from Sohier Dane - Thread Here\n\n

@vaillant identified an issue with the metric that I missed in my initial patch. Essentially, I failed to fully propagate all of the necessary changes when migrating from a pre-launch version of the dataset where the submission and solution file used a wide rather than long format.

\n

The update to the metric is quite simple and is already live. I will be rescoring all existing submissions shortly. That process should take less than an hour.

\n

The metric shouldn't have required multiple patches. I know the bugs have been disruptive and apologize for allowing them to get deployed.

\n

Edit: The rescore is now complete.

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

3  IMPORTS    ⤒

\n\n
\n","metadata":{}},{"cell_type":"code","source":"print(\"\\n... PIP INSTALLS STARTING ...\\n\")\n!pip -qq install pandarallel\nfrom pandarallel import pandarallel\npandarallel.initialize(progress_bar=True)\nprint(\"\\n... PIP INSTALLS COMPLETE ...\\n\")\n\nprint(\"\\n... IMPORTS STARTING ...\\n\")\nprint(\"\\n\\tVERSION INFORMATION\")\n\n# Competition Specific Import\nimport pydicom \nfrom pydicom import dcmread\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut, apply_color_lut, apply_modality_lut\n\nimport pandas as pd; pd.options.mode.chained_assignment = None; pd.set_option('display.max_columns', None);\nimport numpy as np; print(f\"\\t\\t– NUMPY VERSION: {np.__version__}\");\nimport sklearn; print(f\"\\t\\t– SKLEARN VERSION: {sklearn.__version__}\");\nimport cv2; print(f\"\\t\\t– CV2 VERSION: {cv2.__version__}\");\n\n# Built-In Imports (mostly don't worry about these)\nfrom typing import Iterable, Any, Callable, Generator\nfrom kaggle_datasets import KaggleDatasets\nfrom dataclasses import dataclass\nfrom collections import Counter\nfrom datetime import datetime\nfrom zipfile import ZipFile\nfrom glob import glob\nimport subprocess\nimport warnings\nimport requests\nimport textwrap\nimport hashlib\nimport imageio\nimport IPython\nimport urllib\nimport zipfile\nimport pickle\nimport random\nimport shutil\nimport string\nimport json\nimport copy\nimport math\nimport time\nimport gzip\nimport ast\nimport sys\nimport io\nimport gc\nimport re\nimport os\n\n# Visualization Imports (overkill)\nfrom IPython.core.display import HTML, Markdown\nimport matplotlib.pyplot as plt\nfrom matplotlib import animation, rc; rc('animation', html='jshtml')\nfrom tqdm.notebook import tqdm; tqdm.pandas();\nimport plotly.express as px\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 colorama import Fore, Style, init; init()\nimport plotly\nimport PIL\n\ndef hex_to_rgb(hex_color: str) -> tuple:\n \"\"\"Convert hex color to RGB tuple.\n\n Args:\n hex_color (str): The hex color string, starting with '#'.\n\n Returns:\n tuple: A tuple of RGB values.\n \"\"\"\n hex_color = hex_color.lstrip('#')\n return tuple(int(hex_color[i:i+2], 16) for i in (0, 2, 4))\n\ndef clr_print(text: str, color: str = \"#B9508A\", bold: bool = True) -> None:\n \"\"\"Print the given text with the specified color and bold formatting.\n\n Args:\n text (str): The text to format.\n color (str): The hex color code to apply. Defaults to \"#752F55\".\n bold (bool): Whether to apply bold formatting. Defaults to True.\n \"\"\"\n _text = text.replace('\\n', '
')\n rgb = hex_to_rgb(color)\n color_style = f\"color: rgb({rgb[0]}, {rgb[1]}, {rgb[2]});\"\n bold_style = \"font-weight: bold;\" if bold else \"\"\n style = f\"{color_style} {bold_style}\"\n display(HTML(f\"{_text}\"))\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 \nseed_it_all()\n\n# Define the hex colors\nnb_hex_colors = [\"#231942\", \"#5E548E\", \"#9F86C0\", \"#BE95C4\", \"#E0B1CB\", \"#B9508A\", \"#752F55\"]\n\n# Create a Seaborn color palette\nnb_palette = sns.color_palette(nb_hex_colors)\n\n# Is this notebook being run on the backend for scoring re-submission\nIS_DEBUG = False if os.getenv('KAGGLE_IS_COMPETITION_RERUN') else True\nprint(f\"IS DEBUG: {IS_DEBUG}\")\n\n# Plot the palette\nclr_print(\"\\n... NOTEBOOK COLOUR PALETTE ...\")\nsns.palplot(nb_palette)\nplt.show()\n\nprint(\"\\n\\n... IMPORTS COMPLETE ...\\n\")","metadata":{"execution":{"iopub.status.busy":"2024-06-23T19:01:09.472487Z","iopub.execute_input":"2024-06-23T19:01:09.472821Z","iopub.status.idle":"2024-06-23T19:01:27.463483Z","shell.execute_reply.started":"2024-06-23T19:01:09.472792Z","shell.execute_reply":"2024-06-23T19:01:27.462295Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n\n

4  SETUP & HELPER FUNCTIONS    ⤒

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

4.0 FUNCTIONS FROM OTHER KAGGLERS 🧑

\n

\n\nThe first half of the functions below are ripped from the Anatomy & Image Visualization Overview-RSNA RAIDS notebook found here.\n* I went through and converted the relevant code to functions with better documentation, comments, etc.\n* I have also attempted to parameterize the functions to make them more flexible.\n* I have included the code to run, essentially, the whole notebook, as commented out code below the function definitions.\n\nThe second half of the functions come from Cross-Reference Images in Different MRI Planes\n notebook found here.\n* I went through and converted the relevant code to functions with better documentation, comments, etc.\n* I have also attempted to parameterize the functions to make them more flexible.","metadata":{}},{"cell_type":"code","source":"def list_files(directory: str, extension: str = \"\") -> list:\n \"\"\"List all files in a directory with an optional filter for extensions.\n\n Args:\n directory (str): \n Directory path.\n extension (str, optional): \n File extension filter.\n\n Returns:\n list: List of file names.\n \"\"\"\n return [f for f in os.listdir(directory) if f.endswith(extension)]\n\n\ndef plot_diagnosis_distribution(\n df: pd.DataFrame, \n figsize: tuple[int, int] = (20, 5), \n bar_width: float = 0.8\n) -> None:\n \"\"\"Plot the distribution of diagnoses.\n\n Args:\n df (pd.DataFrame): \n DataFrame containing the data.\n figsize (tuple[int, int], optional): \n Size of the figure. Defaults to (20, 5).\n bar_width (float, optional): \n Width of the bars in the bar plot. Defaults to 0.8.\n \"\"\"\n figure, axis = plt.subplots(1, 3, figsize=figsize)\n for idx, d in enumerate(['foraminal', 'subarticular', 'canal']):\n diagnosis = list(filter(lambda x: x.find(d) > -1, df.columns))\n dff = df[diagnosis]\n with warnings.catch_warnings():\n warnings.simplefilter(action='ignore', category=FutureWarning)\n value_counts = dff.apply(pd.value_counts).fillna(0).T\n value_counts.plot(kind='bar', stacked=True, ax=axis[idx], width=bar_width)\n axis[idx].set_title(f'{d} distribution')\n plt.show()\n\n \ndef create_meta_object(part_1: list, base_path: str) -> dict:\n \"\"\"Create a metadata object for the dataset.\n\n Args:\n part_1 (list): \n List of patient directories.\n base_path (str): \n Base path for patient directories.\n\n Returns:\n dict: Metadata object.\n \"\"\"\n meta_obj = {\n p: {\n 'folder_path': os.path.join(base_path, p),\n 'SeriesInstanceUIDs': []\n } for p in part_1\n }\n\n for m in meta_obj:\n meta_obj[m]['SeriesInstanceUIDs'] = list(\n filter(\n lambda x: x.find('.DS') == -1, \n os.listdir(meta_obj[m]['folder_path'])\n )\n )\n return meta_obj\n\n\ndef update_meta_with_descriptions(meta_obj: dict, df_meta: pd.DataFrame) -> None:\n \"\"\"Update the metadata object with series descriptions.\n\n Args:\n meta_obj (dict): \n Metadata object.\n df_meta (pd.DataFrame): \n DataFrame containing series descriptions.\n \"\"\"\n for k in tqdm(meta_obj):\n for s in meta_obj[k]['SeriesInstanceUIDs']:\n if 'SeriesDescriptions' not in meta_obj[k]:\n meta_obj[k]['SeriesDescriptions'] = []\n try:\n series_desc = df_meta[(df_meta['study_id'] == int(k)) & \n (df_meta['series_id'] == int(s))]['series_description'].values[0]\n meta_obj[k]['SeriesDescriptions'].append(series_desc)\n except IndexError:\n meta_obj[k]['SeriesDescriptions'].append('Unknown')\n\n\ndef display_images(\n images: list[np.ndarray], \n title: str, \n cmap: str = 'gray', \n figsize: tuple[int, int] | None = None,\n max_images_per_row: int = 4\n) -> None:\n \"\"\"Display a list of images with a title.\n\n Args:\n images (list[np.ndarray]): \n List of image arrays.\n title (str): \n Title for the plot.\n cmap (str, optional): \n Colormap to use for displaying images. Defaults to 'gray'.\n figsize (tuple[int, int], optional): \n Size of the figure. Defaults to a dynamically calculated size.\n max_images_per_row (int, optional): \n Maximum number of images per row. Defaults to 4.\n \"\"\" \n num_images = len(images)\n num_rows = (num_images + max_images_per_row - 1) // max_images_per_row\n \n if figsize is None:\n figsize= (18, 5*num_rows)\n fig, axes = plt.subplots(num_rows, max_images_per_row, figsize=figsize)\n\n # Flatten axes array if we have multiple rows\n if num_rows > 1:\n axes = axes.flatten()\n else:\n axes = [axes]\n\n for idx, img in enumerate(images):\n ax = axes[idx]\n ax.imshow(img, cmap=cmap)\n ax.axis('off')\n\n # Turn off unused subplots\n for idx in range(num_images, len(axes)):\n axes[idx].axis('off')\n\n plt.suptitle(title)\n plt.tight_layout()\n plt.show()\n\n\ndef display_coor_on_img(\n c: pd.Series, \n img: np.ndarray, \n title: str, \n radius: int = 10, \n color: tuple[int, int, int] = (255, 0, 0), \n thickness: int = 2\n) -> None:\n \"\"\"\n Display coordinates on an image.\n\n Args:\n c (pd.Series): \n Coordinates data.\n img (np.ndarray): \n Image array.\n title (str): \n Title for the plot.\n radius (int, optional): \n Radius of the circle. Defaults to 10.\n color (tuple[int, int, int], optional): \n Color of the circle in BGR format. Defaults to (255, 0, 0).\n thickness (int, optional): \n Thickness of the circle outline. Defaults to 2.\n \"\"\"\n center_coordinates = (int(c['x']), int(c['y']))\n IMG_normalized = cv2.normalize(img, None, alpha=0, beta=255, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_8U)\n IMG_with_circle = cv2.circle(IMG_normalized.copy(), center_coordinates, radius, color, thickness)\n IMG_with_circle = cv2.cvtColor(IMG_with_circle, cv2.COLOR_BGR2RGB)\n plt.imshow(IMG_with_circle)\n plt.axis('off')\n plt.title(title)\n plt.show()\n \n \ndef get_patient_obj(train_df: pd.DataFrame, meta_obj: dict, patient_index: int) -> dict:\n \"\"\"Get the metadata object for a specific patient.\n\n Args:\n train_df (pd.DataFrame): \n DataFrame containing the training data.\n meta_obj (dict): \n Metadata object containing series information.\n patient_index (int): \n Index of the patient in the DataFrame.\n\n Returns:\n dict: \n Metadata object for the specified patient.\n \"\"\"\n patient = train_df.iloc[patient_index]\n patient_obj = meta_obj[str(patient['study_id'])]\n return patient_obj\n\n\ndef get_dicom_image_list(patient_obj: dict) -> dict:\n \"\"\"Get a dictionary of DICOM images for a patient.\n\n Args:\n patient_obj (dict): \n Metadata object for the patient.\n\n Returns:\n dict: \n Dictionary containing DICOM images and their descriptions.\n \"\"\"\n im_list_dcm = {}\n for idx, series_uid in enumerate(patient_obj['SeriesInstanceUIDs']):\n im_list_dcm[series_uid] = {'images': [], 'description': patient_obj['SeriesDescriptions'][idx]}\n images = glob(f\"{patient_obj['folder_path']}/{series_uid}/*.dcm\")\n for img_path in sorted(images, key=lambda x: int(os.path.basename(x).replace('.dcm', ''))):\n im_list_dcm[series_uid]['images'].append({\n 'SOPInstanceUID': os.path.basename(img_path).replace('.dcm', ''), \n 'dicom': pydicom.dcmread(img_path)\n })\n return im_list_dcm\n\n\n# # Load data\n# train_csv_path = '/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train.csv'\n# train = pd.read_csv(train_csv_path)\n# print(f\"Total Cases: {len(train)}\")\n# print(train.columns)\n\n# # Plot diagnosis distribution\n# plot_diagnosis_distribution(train)\n\n# # Process metadata\n# base_path = '/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_images'\n# part_1 = list_files(base_path)\n# meta_obj = create_meta_object(part_1, base_path)\n\n# # Load series descriptions\n# series_desc_csv_path = '/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_series_descriptions.csv'\n# df_meta_f = pd.read_csv(series_desc_csv_path)\n# update_meta_with_descriptions(meta_obj, df_meta_f)\n\n# # Load coordinate data\n# coor_csv_path = '/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_label_coordinates.csv'\n# df_coor = pd.read_csv(coor_csv_path)\n\n# # Get patient object\n# PATIENT_INDEX = 1 # Index of the patient in the training DataFrame\n# patient_obj = get_patient_obj(train, meta_obj, patient_index=PATIENT_INDEX) # Get the metadata object for the specified patient\n# print(patient_obj)\n\n# # Get a dictionary of DICOM images for the patient\n# im_list_dcm = get_dicom_image_list(patient_obj) \n\n# # Display images and coordinates by iterating over each series UID in the DICOM image list\n# for series_uid in im_list_dcm:\n# display_images(\n# [x['dicom'].pixel_array for x in im_list_dcm[series_uid]['images']], # List of DICOM image arrays\n# im_list_dcm[series_uid]['description']\n# )\n\n\n# patient = train.iloc[PATIENT_INDEX] # Get the patient data from the training DataFrame\n# coor_entries = df_coor[df_coor['study_id'] == int(patient['study_id'])] # Filter coordinate entries for the patient\n# print(\"Only showing severe cases for this patient\")\n\n# # Iterate over each coordinate entry and each image in the series\n# for idc, c in coor_entries.iterrows(): \n# for i in im_list_dcm[str(c['series_id'])]['images']: \n# if int(i['SOPInstanceUID']) == int(c['instance_number']): \n \n# # Get the severity level for the patient's condition or set to unknown\n# try: \n# patient_severity = patient[f\"{c['condition'].lower().replace(' ', '_')}_{c['level'].lower().replace('/', '_')}\"] \n# except KeyError:\n# patient_severity = \"unknown severity\" \n \n# # Create the title for the plot\n# title = f\"SOPInstanceUID={i['SOPInstanceUID']}\\nLevel={c['level']}\\nCondition={c['condition']}\\n*SEVERITY={patient_severity}*\\nx={c['x']:.3f}, y={c['y']:.3f}\"\n \n# if patient_severity == 'Severe':\n# # Display the coordinates on the image only if 'Severe'\n# display_coor_on_img(c, i['dicom'].pixel_array, title) \n\n","metadata":{"execution":{"iopub.status.busy":"2024-06-23T19:01:27.465944Z","iopub.execute_input":"2024-06-23T19:01:27.466876Z","iopub.status.idle":"2024-06-23T19:01:27.514040Z","shell.execute_reply.started":"2024-06-23T19:01:27.466836Z","shell.execute_reply":"2024-06-23T19:01:27.512883Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"convert_to_8bit(image: np.ndarray) -> np.ndarray
\nThis function converts a 16-bit image array to an 8-bit image array. This conversion is necessary because many display and image processing tools work best with 8-bit images. The function:\n1. Clips the image values to the 1st and 99th percentiles to reduce the effect of outliers.\n2. Normalizes the image to the range [0, 1].\n3. Scales the normalized values to the range [0, 255] and converts them to `uint8` type.\n\n
\n\nload_dicom_stack(dicom_folder: str, plane: str, reverse_sort: bool = False) -> Dict[str, np.ndarray]
\nThis function loads a stack of DICOM images from a specified folder, sorts them based on their positions, and converts them to 8-bit images.
It returns a dictionary containing:\n1. `array`: The 8-bit image array.\n2. `positions`: The spatial positions of the images.\n3. `pixel_spacing`: The spacing between pixels in the images.\n\nThe `plane` parameter specifies the imaging plane ('sagittal', 'coronal', 'axial').
\nThe `reverse_sort` parameter is used to ensure that the axial images are ordered from caudal to cranial when necessary.\n\n
\n\nplot_dicom_images(study: pd.DataFrame, image_dir: str) -> None
\nThis function plots the middle slices of the sagittal and axial DICOM image stacks from a given study.
\nIt reads the study information from a DataFrame and loads the corresponding DICOM images using the `load_dicom_stack` function.
\nThe images are then displayed side by side for visual comparison.\n\n
\n\nmap_axial_to_sagittal(sag_t2: Dict[str, np.ndarray], ax_t2: Dict[str, np.ndarray]) -> None
\nThis function maps the axial slices to the corresponding sagittal slices and plots the results.
\nIt achieves the following:\n1. Calculates the world space z-coordinates for each y-coordinate in the sagittal image.\n2. Maps the z-coordinates from the axial images to the y-coordinates in the sagittal image.\n3. Plots the sagittal midline slice with red lines indicating the mapped axial slices.\n4. Displays each mapped axial slice alongside the sagittal midline slice for visual inspection.\n\n
\n\n---\n\nCREDIT: The original notebook provided the base code and logic for loading and processing DICOM images, as well as the visualization steps.,br>
This refactored version improves the modularity, readability, and extensibility of the original code by organizing it into well-defined functions with appropriate docstrings and type hinting.\n\n---","metadata":{}},{"cell_type":"code","source":"def convert_to_8bit(image: np.ndarray) -> np.ndarray:\n \"\"\"Convert a 16-bit image to 8-bit.\n \n Args:\n image: \n A 16-bit image array.\n \n Returns:\n An 8-bit image array.\n \"\"\"\n lower, upper = np.percentile(image, (1, 99))\n image = np.clip(image, lower, upper)\n image = image - np.min(image)\n image = image / np.max(image)\n return (image * 255).astype(\"uint8\")\n\n\ndef load_dicom_stack(dicom_folder: str, plane: str, reverse_sort: bool = False) -> dict[str, np.ndarray]:\n \"\"\"Load a stack of DICOM images and convert them to 8-bit.\n \n Args:\n dicom_folder (str): \n Path to the folder containing DICOM files.\n plane (str): \n The plane of the images ('sagittal', 'coronal', 'axial').\n reverse_sort (bool, optional): \n Whether to reverse sort the images.\n \n Returns:\n A dictionary containing the image array, positions, and pixel spacing.\n \"\"\"\n dicom_files = glob(os.path.join(dicom_folder, \"*.dcm\"))\n dicoms = [pydicom.dcmread(f) for f in dicom_files]\n plane_index = {\"sagittal\": 0, \"coronal\": 1, \"axial\": 2}[plane.lower()]\n positions = np.asarray([float(d.ImagePositionPatient[plane_index]) for d in dicoms])\n idx = np.argsort(-positions if reverse_sort else positions)\n ipp = np.asarray([d.ImagePositionPatient for d in dicoms]).astype(\"float\")[idx]\n array = np.stack([d.pixel_array.astype(\"float32\") for d in dicoms])\n array = array[idx]\n return {\n \"array\": convert_to_8bit(array),\n \"positions\": ipp,\n \"pixel_spacing\": np.asarray(dicoms[0].PixelSpacing).astype(\"float\")\n }\n\n\ndef plot_dicom_images(study: pd.DataFrame, image_dir: str) -> None:\n \"\"\"Plot DICOM images for a given study.\n \n Args:\n study (pd.DataFrame): \n DataFrame containing study information.\n image_dir (str): \n Directory containing DICOM images.\n \"\"\"\n for row in study.itertuples():\n if row.series_description == \"Sagittal T2/STIR\":\n sag_t2 = load_dicom_stack(os.path.join(image_dir, str(row.study_id), str(row.series_id)), plane=\"sagittal\")\n elif row.series_description == \"Sagittal T1\":\n sag_t1 = load_dicom_stack(os.path.join(image_dir, str(row.study_id), str(row.series_id)), plane=\"sagittal\")\n elif row.series_description == \"Axial T2\":\n ax_t2 = load_dicom_stack(os.path.join(image_dir, str(row.study_id), str(row.series_id)), plane=\"axial\", reverse_sort=True)\n \n plt.figure(figsize=(14, 5))\n plt.subplot(1, 3, 1)\n plt.imshow(sag_t2[\"array\"][len(sag_t2[\"array\"]) // 2], cmap=\"gray\")\n plt.subplot(1, 3, 2)\n plt.imshow(sag_t1[\"array\"][len(sag_t1[\"array\"]) // 2], cmap=\"gray\")\n plt.subplot(1, 3, 3)\n plt.imshow(ax_t2[\"array\"][len(ax_t2[\"array\"]) // 2], cmap=\"gray\")\n plt.show()\n\n\ndef map_axial_to_sagittal(sag_t2: dict[str, np.ndarray], ax_t2: dict[str, np.ndarray], plot_all_slices: bool = True) -> None:\n \"\"\"Map axial slices to sagittal slices and plot the results.\n \n Args:\n sag_t2 (dict[str, np.ndarray]): \n Dictionary containing sagittal T2 image data.\n ax_t2 (dict[str, np.ndarray]): \n Dictionary containing axial T2 image data.\n plot_all_slices (bool, optional):\n Whether to show each individual side-by-side plot.\n \"\"\"\n top_left_hand_corner_sag_t2 = sag_t2[\"positions\"][len(sag_t2[\"array\"]) // 2]\n sag_y_axis_to_pixel_space = [top_left_hand_corner_sag_t2[2]]\n while len(sag_y_axis_to_pixel_space) < sag_t2[\"array\"].shape[1]:\n sag_y_axis_to_pixel_space.append(sag_y_axis_to_pixel_space[-1] - sag_t2[\"pixel_spacing\"][1])\n \n sag_y_coord_to_axial_slice = {}\n for ax_t2_slice, ax_t2_pos in zip(ax_t2[\"array\"], ax_t2[\"positions\"]):\n diffs = np.abs(np.asarray(sag_y_axis_to_pixel_space) - ax_t2_pos[2])\n sag_y_coord = np.argmin(diffs)\n sag_y_coord_to_axial_slice[sag_y_coord] = ax_t2_slice\n \n sag_midline_slice = sag_t2[\"array\"][len(sag_t2[\"array\"]) // 2]\n plt.imshow(sag_midline_slice, cmap=\"gray\")\n for k in [*sag_y_coord_to_axial_slice]:\n plt.axhline(y=k, color=\"red\", linestyle=\"--\")\n plt.show()\n \n if plot_all_slices:\n for k, v in sag_y_coord_to_axial_slice.items():\n plt.figure(figsize=(10, 4))\n plt.subplot(1, 2, 1)\n plt.imshow(sag_midline_slice, cmap=\"gray\")\n plt.axhline(y=k, color=\"red\", linestyle=\"--\")\n plt.subplot(1, 2, 2)\n plt.imshow(v, cmap=\"gray\")\n plt.show()\n \n\n \n_df = pd.read_csv(\"/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_series_descriptions.csv\")\n_image_dir = \"/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_images/\"\nstudy = _df.loc[_df.study_id == _df.study_id.iloc[10]]\nplot_dicom_images(study, _image_dir)\nax_t2 = sag_t1 = sag_t2 = None \nfor row in study.itertuples():\n if row.series_description == \"Sagittal T2/STIR\":\n sag_t2 = load_dicom_stack(os.path.join(_image_dir, str(row.study_id), str(row.series_id)), plane=\"sagittal\")\n elif row.series_description == \"Sagittal T1\":\n sag_t1 = load_dicom_stack(os.path.join(_image_dir, str(row.study_id), str(row.series_id)), plane=\"sagittal\")\n elif row.series_description == \"Axial T2\":\n ax_t2 = load_dicom_stack(os.path.join(_image_dir, str(row.study_id), str(row.series_id)), plane=\"axial\", reverse_sort=True)\n \n if sag_t2 and ax_t2:\n break\n\n# Assuming the variables sag_t2 and ax_t2 are defined in the study\nmap_axial_to_sagittal(sag_t2, ax_t2)\n\ndel _df, _image_dir, sag_t1, sag_t2, ax_t2, study; gc.collect(); gc.collect();","metadata":{"execution":{"iopub.status.busy":"2024-06-23T19:01:27.520570Z","iopub.execute_input":"2024-06-23T19:01:27.524596Z","iopub.status.idle":"2024-06-23T19:01:39.024654Z","shell.execute_reply.started":"2024-06-23T19:01:27.524540Z","shell.execute_reply":"2024-06-23T19:01:39.023461Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

4.1 GENERIC FUNCTIONS

\n

\n\nSome functions I bring with me...","metadata":{}},{"cell_type":"code","source":"def flatten_l_o_l(nested_list):\n \"\"\" Flatten a list of lists into a single list.\n\n Args:\n nested_list (Iterable): \n – A list of lists (or iterables) to be flattened.\n\n Returns:\n A flattened list containing all items from the input list of lists.\n \"\"\"\n return [item for sublist in nested_list for item in sublist]\n\n\ndef print_ln(symbol=\"-\", line_len=110, newline_before=False, newline_after=False):\n \"\"\" Print a horizontal line of a specified length and symbol.\n\n Args:\n symbol (str, optional): \n – The symbol to use for the horizontal line\n line_len (int, optional): \n – The length of the horizontal line in characters\n newline_before (bool, optional): \n – Whether to print a newline character before the line\n newline_after (bool, optional): \n – Whether to print a newline character after the line\n \n Returns:\n None; A divider with pre/post new-lines (optional) is printed\n \"\"\"\n if newline_before: print();\n print(symbol * line_len)\n if newline_after: print();\n \n \ndef display_hr(newline_before=False, newline_after=False):\n \"\"\" Renders a HTML
\n\n Args:\n newline_before (bool, optional): \n – Whether to print a newline character before the line\n newline_after (bool, optional): \n – Whether to print a newline character after the line\n \n Returns:\n None; A divider with pre/post new-lines (optional) is printed\n \"\"\"\n if newline_before: print();\n display(HTML(\"
\"))\n if newline_after: print();\n\n\ndef wrap_text(text, width=88):\n \"\"\"Wrap text to a specified width.\n\n Args:\n text (str): \n - The text to wrap.\n width (int): \n - The maximum width of a line. Default is 88.\n\n Returns:\n str: The wrapped text.\n \"\"\"\n return textwrap.fill(text, width)\n\n\ndef wrap_text_by_paragraphs(text, width=88):\n \"\"\"Wrap text by paragraphs to a specified width.\n\n Args:\n text (str): \n - The text containing multiple paragraphs to wrap.\n width (int): \n - The maximum width of a line. Default is 88.\n\n Returns:\n str: The wrapped text with preserved paragraph separation.\n \"\"\"\n paragraphs = text.split('\\n') # Assuming paragraphs are separated by newlines\n wrapped_paragraphs = [textwrap.fill(paragraph, width) for paragraph in paragraphs]\n return '\\n\\n'.join(wrapped_paragraphs)","metadata":{"execution":{"iopub.status.busy":"2024-06-23T19:01:39.027293Z","iopub.execute_input":"2024-06-23T19:01:39.028261Z","iopub.status.idle":"2024-06-23T19:01:39.041478Z","shell.execute_reply.started":"2024-06-23T19:01:39.028221Z","shell.execute_reply":"2024-06-23T19:01:39.039776Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

4.2 COMPETITION HELPERS

\n

\n\n","metadata":{}},{"cell_type":"code","source":"def dicom_path_to_image(path, try_lut=True, fix_monochrome=True):\n \"\"\" Convert dicom file to numpy array \n \n Args:\n path (str): \n Path to the dicom file to be converted\n try_lut (bool): \n Whether or not VOI LUT is available.\n VOI LUT (if available by DICOM device) is used to transform raw DICOM data to \"human-friendly\" view\n fix_monochrome (bool): \n Whether or not to apply monochrome fix\n \n Returns:\n Numpy array of the respective dicom file \n \n \"\"\"\n # (1) Use the pydicom library to read the dicom file and get image array\n dicom = pydicom.read_file(path)\n arr = dicom.pixel_array\n \n # (2) Some DICOM datasets store their output image pixel values in a lookup table (LUT)\n # - The values in Pixel Data are the index to a corresponding LUT entry. \n # - When a dataset’s (0028,0004) Photometric Interpretation value is PALETTE COLOR then we should ...\n # use the apply_color_lut() function to apply a palette color LUT to the pixel data to produce an RGB image.\n if try_lut and dicom.PhotometricInterpretation==\"PALETTE COLOR\":\n clr_print(\"\\n\\n... Applying COLOR LUT ...\\n\\n\")\n arr = apply_color_lut(arr, ds)\n \n \n # (3) The DICOM Modality LUT module (similar to the Color one) converts raw pixel data values to a specific (possibly unitless) physical quantity.\n # Examples are quantities such as Hounsfield units for CT scan . \n # The apply_modality_lut() function can be used with an input array of raw values and a dataset containing a Modality LUT module to return the converted values. \n # When a dicom dataset requires multiple grayscale transformations, the Modality LUT transformation is always applied first.\n if try_lut:\n hu = apply_modality_lut(arr, dicom)\n if (arr!=hu).any(): \n clr_print(\"\\n\\n... Applying MODALITY LUT ...\\n\\n\")\n \n # (4) The DICOM VOI LUT module applies a VOI or windowing operation to input values. \n # The apply_voi_lut() function can be used with an input array and a dataset containing a VOI LUT module to return values with applied VOI LUT or windowing. \n # When a dicom dataset contains multiple VOI or windowing views then a particular view can be returned by using the index keyword parameter. \n # In this case the index 0 will be used.\n # When a dataset requires multiple grayscale transformations, then it’s assumed that the modality LUT or rescale operation has already been applied.\n if try_lut:\n arr = apply_voi_lut(hu, dicom, index=0)\n if (arr!=hu).any(): clr_print(\"\\n\\n... Applying VOI LUT ...\\n\\n\")\n \n # The XRAY may look inverted\n # - If we want to fix this we can\n if fix_monochrome and dicom.PhotometricInterpretation == \"MONOCHROME1\":\n clr_print(\"\\n\\n... Applying MONOCHROME FIX ...\\n\\n\")\n arr = np.amax(arr) - arr\n \n # Normalize the image array and return\n lower, upper = np.percentile(x, (1, 99))\n arr = np.clip(arr, lower, upper)\n arr = arr - np.min(arr)\n arr = arr / np.max(arr)\n arr = (arr * 255).astype(np.uint8)\n return arr\n\n\ndef dicom_array_to_image(dicom, arr, try_lut=True, fix_monochrome=True):\n \"\"\"Convert dicom array to numpy image array.\n \n Args:\n dicom (pydicom.dataset.FileDataset): \n DICOM object containing metadata\n dicom_array (numpy.ndarray): \n Numpy array containing DICOM pixel data\n try_lut (bool): \n Whether or not VOI LUT is available.\n VOI LUT (if available by DICOM device) is used to transform raw DICOM data to \"human-friendly\" view\n fix_monochrome (bool): \n Whether or not to apply monochrome fix\n \n Returns:\n Numpy array of the respective dicom file \n \"\"\"\n if try_lut and dicom.PhotometricInterpretation == \"PALETTE COLOR\":\n arr = apply_color_lut(arr, dicom)\n \n if try_lut:\n arr = apply_modality_lut(arr, dicom)\n \n if try_lut:\n arr = apply_voi_lut(arr, dicom, index=0)\n \n if fix_monochrome and dicom.PhotometricInterpretation == \"MONOCHROME1\":\n arr = np.amax(arr) - arr\n \n lower, upper = np.percentile(arr, (1, 99))\n arr = np.clip(arr, lower, upper)\n arr = arr - np.min(arr)\n arr = arr / np.max(arr)\n arr = (arr * 255).astype(np.uint8)\n return arr","metadata":{"execution":{"iopub.status.busy":"2024-06-23T19:01:39.043256Z","iopub.execute_input":"2024-06-23T19:01:39.044094Z","iopub.status.idle":"2024-06-23T19:01:39.064151Z","shell.execute_reply.started":"2024-06-23T19:01:39.044028Z","shell.execute_reply":"2024-06-23T19:01:39.062732Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"_POSSIBLE_DICOM_ATTRS = [\n \"BitsAllocated\", \"BitsStored\", \"Columns\", \"ContentDate\", \"ContentTime\", \"FrameOfReferenceUID\", \n \"HighBit\", \"ImageOrientationPatient\", \"ImagePositionPatient\", \"InstanceNumber\", \n \"PatientID\", \"PatientPosition\", \"PhotometricInterpretation\", \"PixelData\", \n \"PixelRepresentation\", \"PixelSpacing\", \"RescaleIntercept\", \"RescaleSlope\", \"RescaleType\", \n \"Rows\", \"SOPInstanceUID\", \"SamplesPerPixel\", \"SeriesDescription\", \"SeriesInstanceUID\", \n \"SliceLocation\", \"SliceThickness\", \"SpacingBetweenSlices\", \"StudyInstanceUID\", \n \"WindowCenter\", \"WindowWidth\"\n]\n\ndef extract_all_dcm_data(dcm_path: str, save_to_dir: str = \"/kaggle/working/pngs/train\", save_to_png: bool = True) -> dict:\n \"\"\"Extract all DICOM data and optionally save the image to PNG.\n \n Args:\n dcm_path (str): Path to the DICOM file.\n save_to_dir (str): Directory to save the PNG images.\n save_to_png (bool): Whether or not to save the image as PNG.\n \n Returns:\n Dicom attributes\n \"\"\"\n dicom = pydicom.read_file(dcm_path) \n dicom_attr_dict = {\n attr_key:dicom.get(attr_key)\n for attr_key in _POSSIBLE_DICOM_ATTRS\n if attr_key!=\"PixelData\"\n }\n \n if save_to_png:\n # Get save path info and create if not existing\n file_path_ending = \"/\".join(dcm_path.rsplit(\"/\", 3)[1:]).replace(\".dcm\", \".png\")\n save_path = os.path.join(save_to_dir, file_path_ending)\n os.makedirs(save_path.rsplit(\"/\", 1)[0], exist_ok=True)\n \n # Save the image\n img = Image.fromarray(dicom_array_to_image(dicom, dicom.pixel_array))\n img.save(save_path)\n dicom_attr_dict[\"PNGPath\"] = save_path\n \n return dicom_attr_dict\n\n\ndef create_dicom_df(all_dcm_paths: list[str], is_train: bool = True, save_to_png: bool = True) -> pd.DataFrame:\n \"\"\"Create a DataFrame with DICOM data and optionally save images as PNGs.\n \n Args:\n all_dcm_paths (list[str]): List of paths to DICOM files.\n is_train (bool): Whether the data is training data or not.\n save_to_dir (str): Directory to save the PNG images.\n save_to_png (bool): Whether or not to save the image as PNG.\n \n Returns:\n pd.DataFrame: DataFrame containing DICOM attributes and paths.\n \"\"\"\n # Set PNG directory based on training or test data\n _png_dir = \"/kaggle/working/pngs/train\" if is_train else \"/kaggle/working/pngs/test\"\n \n # Create the initial DataFrame with DICOM paths\n dicom_df = pd.DataFrame({\"dcm_path\": all_dcm_paths})\n \n # Extract study_id, series_id, and instance_number from DICOM paths\n dicom_df[[\"study_id\", \"series_id\", \"instance_number\"]] = pd.DataFrame(\n dicom_df.dcm_path.apply(lambda x: [x.replace(\".dcm\", \"\") for x in x.rsplit(\"/\", 3)[1:]]).tolist()\n ).astype(\"int\")\n \n # Extract DICOM attributes and optionally save images as PNGs\n _new_df_cols = [x for x in _POSSIBLE_DICOM_ATTRS if x != \"PixelData\"]\n if save_to_png:\n _new_df_cols += [\"PNGPath\",]\n \n # Create the new dataframe using parallel processing...\n # I did this on Paperspace with 32 CPUs (vs the 2 on Kaggle)\n dicom_df[_new_df_cols] = pd.DataFrame(dicom_df[\"dcm_path\"].parallel_apply(lambda x: extract_all_dcm_data(x, save_to_dir=_png_dir)).tolist())\n \n return dicom_df\n\ndef get_study_labels(df: pd.DataFrame, study_id: int | str, one_hot_encode_labels: bool = False, sparse_encode_labels: bool = False,) -> dict[str, str | int | list[float]]:\n \"\"\"Get the labels for a given study from a pandas dataframe (train).\n\n Args:\n df (pd.DataFrame): \n DataFrame containing study data.\n study_id (int | str): \n ID of the study to retrieve labels for.\n one_hot_encode_labels (bool, optional): \n Whether to one-hot encode string labels.\n sparse_encode_labels (bool, optional): \n Whether to encode string labels as integers.\n \n Raises:\n ValueError: \n If study_id is not found within the provided dataframe.\n \n Returns:\n dict[str, str | int | list[float]]: \n Dictionary with columns as keys and labels (either strings, integers, or one-hot encoded lists) as values.\n \"\"\"\n # Extract the row corresponding to the given study_id\n study_row = df[df['study_id'] == study_id]\n \n # Ensure the study_id exists in the dataframe\n if study_row.empty:\n raise ValueError(f\"Study ID {study_id} not found in the dataframe.\")\n\n # Convert the single-row DataFrame to a dictionary\n labels_dict = study_row.iloc[0].to_dict()\n\n # Remove 'study_id' from the dictionary\n labels_dict.pop('study_id')\n\n # Encode labels if requested (priority given to one_hot)\n if one_hot_encode_labels:\n labels_dict = {key: [1.0 if i == str2int_severity[value] else 0.0 for i in range(len(LABEL_STRS))] for key, value in labels_dict.items()}\n elif sparse_encode_labels:\n labels_dict = {key: str2int_severity[value] for key, value in labels_dict.items()}\n return labels_dict","metadata":{"execution":{"iopub.status.busy":"2024-06-23T19:01:39.066476Z","iopub.execute_input":"2024-06-23T19:01:39.066942Z","iopub.status.idle":"2024-06-23T19:01:39.087739Z","shell.execute_reply.started":"2024-06-23T19:01:39.066893Z","shell.execute_reply":"2024-06-23T19:01:39.086262Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

4.3 DATASET LOADING

\n

\n\nLet's load the dataset our own way and improve upon the previously shown code","metadata":{}},{"cell_type":"code","source":"# ROOT PATHS\nWORKING_DIR = \"/kaggle/working\"\nINPUT_DIR = \"/kaggle/input\"\nCOMPETITION_DIR = os.path.join(INPUT_DIR, \"rsna-2024-lumbar-spine-degenerative-classification\")\nTRAIN_DCM_IMAGE_DIR = os.path.join(COMPETITION_DIR, \"train_images\")\nTEST_DCM_IMAGE_DIR = os.path.join(COMPETITION_DIR, \"test_images\")\nADDITIONAL_INPUT_DIR = \"/kaggle/input/rsna-lsdc-files/rsna_lsdc_files\"\nTRAIN_DICOM_ATTR_CSV = os.path.join(ADDITIONAL_INPUT_DIR, \"train_dicom.csv\")\nTEST_DICOM_ATTR_CSV = os.path.join(ADDITIONAL_INPUT_DIR, \"train_dicom.csv\")\nTRAIN_PNG_IMAGE_DIR = os.path.join(ADDITIONAL_INPUT_DIR, \"pngs\", \"train\")\nTRAIN_PNG_IMAGE_DIR = os.path.join(ADDITIONAL_INPUT_DIR, \"pngs\", \"test\")\n\n\n# COMPETITION FILE PATHS\nSS_CSV_PATH = os.path.join(COMPETITION_DIR, \"sample_submission.csv\")\nTRAIN_CSV_PATH = os.path.join(COMPETITION_DIR, \"train.csv\")\nTRAIN_SERIES_DESC_CSV_PATH = os.path.join(COMPETITION_DIR, \"train_series_descriptions.csv\")\nTRAIN_LABEL_COORDINATES_CSV_PATH = os.path.join(COMPETITION_DIR, \"train_label_coordinates.csv\")\nTEST_SERIES_DESC_CSV_PATH = os.path.join(COMPETITION_DIR, \"test_series_descriptions.csv\")\n\n# DEFINE COMPETITION DATAFRAMES\nclr_print(\"\\n\\n... SAMPLE SUBMISSION DATAFRAME ...\\n\\n\")\nss_df = pd.read_csv(SS_CSV_PATH)\ndisplay(ss_df)\n\nclr_print(\"\\n\\n... TRAIN DATAFRAME ...\\n\\n\")\ntrain_df = pd.read_csv(TRAIN_CSV_PATH)\ndisplay(train_df)\n\nclr_print(\"\\n\\n... TRAIN SERIES DESCRIPTIONS ...\\n\\n\")\ntrain_series_desc_df = pd.read_csv(TRAIN_SERIES_DESC_CSV_PATH)\ndisplay(train_series_desc_df)\n\nclr_print(\"\\n\\n... TEST SERIES DESCRIPTIONS ...\\n\\n\")\ntest_series_desc_df = pd.read_csv(TEST_SERIES_DESC_CSV_PATH)\ndisplay(test_series_desc_df)\n\nclr_print(\"\\n\\n... TRAIN LABEL COORDINATES ...\\n\\n\")\ntrain_label_coords_df = pd.read_csv(TRAIN_LABEL_COORDINATES_CSV_PATH)\ndisplay(train_label_coords_df)\n\n# DEFINE CONSTANTS AND USEFUL GLOBALS\nSERIES_TYPES = [\"Axial T2\", \"Sagittal T1\", \"Sagittal T2/STIR\"]\nLABEL_STRS = ['Normal/Mild', 'Moderate', 'Severe']\nstr2int_severity = {lbl:i for i, lbl in enumerate(LABEL_STRS)}\nint2str_severity = {v:k for k,v in str2int_severity.items()}\nLABEL_INTS = [str2int_severity[lbl] for lbl in LABEL_STRS]\n\nDICOM_NAMING_MAP = {\n \"PNGPath\": \"png_path\", \"Columns\":\"img_width\", \"Rows\":\"img_height\", \n \"SliceLocation\": \"slice_location\", \"SliceThickness\": \"slice_thickness\", \n \"PixelSpacing\": \"pixel_spacing\", \"SpacingBetweenSlices\": \"slice_spacing\",\n \"WindowCenter\": \"window_center\", \"WindowWidth\": \"window_width\", \n}\n\n# Redundant, already handled in loading/png-creation, or not enough variability to retain\nDICOM_COLS_TO_DROP = [\n \"StudyInstanceUID\", \"SeriesInstanceUID\", \"InstanceNumber\", # Redundant\n \"BitsAllocated\", \"BitsStored\", \"HighBit\", \"PixelRepresentation\", # These determine how pixel values are stored and their range. The current loading/saving function should handle this already (I THINK!) --> https://pydicom.github.io/pydicom/stable/old/image_data_handlers.html\n \"ContentDate\", \"ContentTime\", # These fields are generally not relevant unless the model specifically requires temporal context.\n \"FrameOfReferenceUID\", \"SOPInstanceUID\", # Typically used for ensuring data integrity and traceability in clinical settings but not directly for image processing (also often 1-to-1 with series_id/instance_id)\n \"PhotometricInterpretation\", # Handled by our loading function,\n \"PatientID\", # Always 1-to-1 with study-id -- so it is superfluous\n \"RescaleIntercept\", # Only 1 value - Always 0.0 when present.\n \"PhotometricInterpretation\", # Always MONOCHROME2\n]\nif not os.path.isfile(TRAIN_DICOM_ATTR_CSV):\n # CREATE A DICOM DATAFRAME... AND OPTIONALLY SAVE THE PIXEL ARRAYS AS PNGS\n ALL_TRAIN_DCM_PATHS = glob(os.path.join(TRAIN_DCM_IMAGE_DIR, \"**\", \"*.dcm\"), recursive=True)\n train_dicom_df = create_dicom_df(ALL_TRAIN_DCM_PATHS, save_to_png=False)\nelse:\n train_dicom_df = pd.read_csv(TRAIN_DICOM_ATTR_CSV)\n train_dicom_df[\"dcm_path\"] = train_dicom_df[\"dcm_path\"].str.replace(WORKING_DIR, ADDITIONAL_INPUT_DIR)\n ALL_TRAIN_DCM_PATHS = train_dicom_df[\"dcm_path\"].tolist()\n train_dicom_df[\"PNGPath\"] = train_dicom_df[\"PNGPath\"].str.replace(WORKING_DIR, ADDITIONAL_INPUT_DIR)\ntrain_dicom_df = pd.merge(train_series_desc_df, train_dicom_df, on=[\"study_id\", \"series_id\"])\ntrain_dicom_df = train_dicom_df.rename(columns=DICOM_NAMING_MAP).drop(columns=DICOM_COLS_TO_DROP)\ntrain_dicom_df = train_dicom_df[[x for x in train_dicom_df.columns if \"_\" in x]+[x for x in train_dicom_df.columns if \"_\" not in x]]\n\nALL_TEST_DCM_PATHS = glob(os.path.join(TEST_DCM_IMAGE_DIR, \"**\", \"*.dcm\"), recursive=True)\ntest_dicom_df = create_dicom_df(ALL_TEST_DCM_PATHS, save_to_png=True)\ntest_dicom_df = pd.merge(test_series_desc_df, test_dicom_df, on=[\"study_id\", \"series_id\"])\ntest_dicom_df = test_dicom_df.rename(columns=DICOM_NAMING_MAP).drop(columns=DICOM_COLS_TO_DROP)\ntest_dicom_df = test_dicom_df[[x for x in test_dicom_df.columns if \"_\" in x]+[x for x in test_dicom_df.columns if \"_\" not in x]]\n\nclr_print(\"\\n\\n... TRAIN DICOM DF ...\\n\\n\")\ndisplay(train_dicom_df)\n\nclr_print(\"\\n\\n... TEST DICOM DF ...\\n\\n\")\ndisplay(test_dicom_df)","metadata":{"execution":{"iopub.status.busy":"2024-06-23T19:01:39.089369Z","iopub.execute_input":"2024-06-23T19:01:39.090083Z","iopub.status.idle":"2024-06-23T19:01:47.219582Z","shell.execute_reply.started":"2024-06-23T19:01:39.090042Z","shell.execute_reply":"2024-06-23T19:01:47.218460Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def simplify_orientation(orientation):\n \"\"\" TBD \"\"\" \n # Convert from string to array of floats\n if isinstance(orientation, str):\n orientation = ast.literal_eval(orientation)\n \n row_cosines, col_cosines = orientation[:3], orientation[3:]\n \n # Calculate the dominant axis\n row_dominant_axis = np.argmax(np.abs(row_cosines))\n col_dominant_axis = np.argmax(np.abs(col_cosines))\n \n # Compute the angle between vectors\n angle = np.arccos(np.dot(row_cosines, col_cosines) / (np.linalg.norm(row_cosines) * np.linalg.norm(col_cosines)))\n \n return row_dominant_axis, col_dominant_axis, np.round(np.degrees(angle), 2)\n\ntrain_dicom_df[\"primary_axes_with_relative_angle\"] = train_dicom_df[\"ImageOrientationPatient\"].progress_apply(simplify_orientation)\ntrain_dicom_df = train_dicom_df.drop(columns=[\"ImageOrientationPatient\"])\n\n####\n# Occurence Are Below... NOTE: all relative angles are ~90 degrees indicating orthogonality.\n####\n# (0, 1, 90.0): \n# - This is the most common orientation\n# - Primary alignment is as follows:\n# - the x-axis (left-right) and \n# - the y-axis (anterior-posterior)\n# - This corresponds to the axial plane. \n# (1, 2, 90.0): \n# - This is the second most common (almost as common as Axial) orientation. \n# - Primary alignment is as follows: \n# - the y-axis (anterior-posterior) and \n# - the z-axis (superior-inferior)\n# - This corresponds to the sagittal plane. \n# (0, 2, 90.0): \n# - Rare and possibly erroneous.\n# - Primary alignment is as follows:\n# - The x-axis (left-right) and \n# - The z-axis (superior-inferior)].\n# - This could either be a variant of sagittal scans or another less common plane like coronal (?)\n####\n\ntrain_dicom_df","metadata":{"execution":{"iopub.status.busy":"2024-06-23T19:01:47.221322Z","iopub.execute_input":"2024-06-23T19:01:47.222232Z","iopub.status.idle":"2024-06-23T19:01:57.692173Z","shell.execute_reply.started":"2024-06-23T19:01:47.222196Z","shell.execute_reply":"2024-06-23T19:01:57.691007Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**attempting to plot orientation and patient location... probably unnecessary**","metadata":{}},{"cell_type":"code","source":"# import plotly.graph_objects as go\n\n# def plot_orientations(orientations: list, positions: list = None, title: str = \"Image Orientations\") -> None:\n# \"\"\"\n# Plot multiple image orientations and their positions in 3D using Plotly.\n\n# Args:\n# orientations (list): List of orientations, each containing six floats \n# representing two direction cosine vectors.\n# positions (list, optional): List of positions, each containing three floats \n# representing the (X, Y, Z) coordinates. Defaults to None.\n# title (str): Title of the plot. Defaults to \"Image Orientations\".\n \n# Example:\n# orientations = [\n# [1.0E-13, 1.0, 4.896E-12, -0.020486, 4.897E-12, -0.99979],\n# [-0, 1, 0, -0, -0, -1],\n# [-3.222e-012, 1, -2.05078e-010, -0.0157073092132, -2.05103e-010, -0.9998766326089]\n# ]\n# positions = [\n# [16.8411, -108.852, 131.582],\n# [-1.7507368326, -117.9880981445, 224.0385894775],\n# [14.842441574236, -66.392250032321, 139.60760050879]\n# ]\n# plot_orientations(orientations, positions)\n# \"\"\"\n# fig = go.Figure()\n\n# # Determine the maximum coordinate value for scaling\n# if positions is None:\n# positions = [[0, 0, 0] for _ in orientations]\n \n# xs, ys, zs = [p[0] for p in positions], [p[1] for p in positions], [p[2] for p in positions]\n# coord_extents = ((min(xs)-1, max(xs)+1), (min(ys)-1, max(ys)+1), (min(zs)-1, max(zs)+1))\n# max_coord_magnitude = abs(max([_max-_min for (_min, _max) in coord_extents]))\n# # Add the patient coordinate system\n# fig.add_trace(go.Scatter3d(x=[coord_extents[0][0], coord_extents[0][1]], y=[0, 0], z=[0, 0], mode='lines', line=dict(color='red', width=5), name='X-axis'))\n# fig.add_trace(go.Scatter3d(x=[0, 0], y=[coord_extents[1][0], coord_extents[1][1]], z=[0, 0], mode='lines', line=dict(color='green', width=5), name='Y-axis'))\n# fig.add_trace(go.Scatter3d(x=[0, 0], y=[0, 0], z=[coord_extents[2][0], coord_extents[2][1]], mode='lines', line=dict(color='blue', width=5), name='Z-axis'))\n\n# for idx, (orientation, position) in enumerate(zip(orientations, positions)):\n# row_direction = orientation[:3]\n# col_direction = orientation[3:]\n\n# # Plot patient position as a dot\n# fig.add_trace(go.Scatter3d(x=[position[0]], y=[position[1]], z=[position[2]], \n# mode='markers', marker=dict(size=5, color='black'), name=f'Position {idx+1}'))\n\n\n# # Calculate the end points of the vectors\n# vec_end = [position[i] + row_direction[i] * max_coord_magnitude + col_direction[i] * max_coord_magnitude for i in range(3)]\n\n# # Plot row and column direction lines\n# fig.add_trace(go.Scatter3d(x=[position[0], vec_end[0]], \n# y=[position[1], vec_end[1]], \n# z=[position[2], vec_end[2]], \n# mode='lines', line=dict(color='cyan', width=5), name=f'R{idx+1}'))\n\n# # Add arrowheads using cones\n# fig.add_trace(go.Cone(x=[vec_end[0]], y=[vec_end[1]], z=[vec_end[2]],\n# u=[row_direction[0]*max_coord_magnitude+col_direction[0]*max_coord_magnitude], \n# v=[row_direction[1]*max_coord_magnitude+col_direction[1]*max_coord_magnitude], \n# w=[row_direction[2]*max_coord_magnitude+col_direction[2]*max_coord_magnitude],\n# colorscale=[[0, 'cyan'], [1, 'cyan']], sizemode='scaled', sizeref=0.1,\n# showscale=False, \n# anchor=\"cm\"\n# )\n# )\n\n# # Setting labels and title\n# fig.update_layout(scene=dict(\n# xaxis=dict(title='X', range=[coord_extents[0][0], coord_extents[0][1]]),\n# yaxis=dict(title='Y', range=[coord_extents[1][0], coord_extents[1][1]]),\n# zaxis=dict(title='Z', range=[coord_extents[2][0], coord_extents[2][1]])),\n# title=title\n# )\n\n# fig.show()\n\n\n\n# plot_orientations(orientations, positions)\n\n# for x, _df in train_dicom_df.groupby(\"series_id\"):\n# if x==10996:\n# continue\n# orientations = _df[\"ImageOrientationPatient\"].apply(ast.literal_eval).tolist()\n# positions = _df[\"ImagePositionPatient\"].apply(ast.literal_eval).tolist()\n \n# plot_orientations(orientations, positions, title=f\"PATIENT {x}\")\n# break","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2024-06-23T19:01:57.693833Z","iopub.execute_input":"2024-06-23T19:01:57.694173Z","iopub.status.idle":"2024-06-23T19:01:57.705248Z","shell.execute_reply.started":"2024-06-23T19:01:57.694143Z","shell.execute_reply":"2024-06-23T19:01:57.703846Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Merging series descriptions with label coordinates\ntrain_df_series_coords = pd.merge(train_series_desc_df, train_label_coords_df, on=['study_id', 'series_id'])\n\n# Standardize\ntrain_df_series_coords['condition'] = train_df_series_coords['condition'].str.replace(\" \", \"_\").str.lower()\ntrain_df_series_coords['level'] = train_df_series_coords['level'].str.upper().str.replace('_', '/')\n\n# Reshape the train dataframe to long format\ntrain_df_long = train_df.melt(id_vars=['study_id'], var_name='condition_level', value_name='severity')\ntrain_df_long['condition'] = train_df_long['condition_level'].apply(lambda x: x.rsplit(\"_\", 2)[0])\ntrain_df_long['level'] = train_df_long['condition_level'].apply(lambda x: \"/\".join(x.rsplit(\"_\", 2)[-2:]).upper().replace(\"_\", \"/\"))\ntrain_df_long = train_df_long.drop('condition_level', axis=1).sort_values(by=[\"study_id\"]).reset_index(drop=True)\n\n# Combine the reshaped train_df_long with train_df_series_coords based on study_id, condition, and level.\ntrain_df_combined = pd.merge(train_df_series_coords, train_df_long, on=['study_id', 'condition', 'level'])\ntrain_df_combined[\"dcm_path\"] = (\n TRAIN_DCM_IMAGE_DIR +\"/\"+\n train_df_combined[\"study_id\"].astype(str) +\"/\"+\n train_df_combined[\"series_id\"].astype(str) +\"/\"+\n train_df_combined[\"instance_number\"].astype(str) +\".dcm\"\n)\n\nclr_print(\"\\n\\n... MELTED/COMBINED TRAINING DATAFRAME (1 ROW PER CONDITION/LEVEL) ...\\n\\n\")\ntrain_df_combined","metadata":{"execution":{"iopub.status.busy":"2024-06-23T19:01:57.708404Z","iopub.execute_input":"2024-06-23T19:01:57.708764Z","iopub.status.idle":"2024-06-23T19:01:58.027582Z","shell.execute_reply.started":"2024-06-23T19:01:57.708734Z","shell.execute_reply":"2024-06-23T19:01:58.026501Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DEMO_STUDY_ID = 4646740\nDEMO_STUDY_LABEL_MAP = get_study_labels(train_df, DEMO_STUDY_ID)\nDEMO_DICOM_DF = train_dicom_df[train_dicom_df.study_id==DEMO_STUDY_ID].reset_index(drop=True)\nDEMO_SERIES_IDS = DEMO_DICOM_DF.series_id.unique().tolist()\n\n# SET THIS UP AS A CLASS\n# def dicom_df__plot_series(df):\nDEMO_DICOM_DF.series_description.value_counts()","metadata":{"execution":{"iopub.status.busy":"2024-06-23T19:01:58.029175Z","iopub.execute_input":"2024-06-23T19:01:58.029638Z","iopub.status.idle":"2024-06-23T19:01:58.044347Z","shell.execute_reply.started":"2024-06-23T19:01:58.029598Z","shell.execute_reply":"2024-06-23T19:01:58.043151Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def dicom_df__get_series_image_path(dicom_df: pd.DataFrame, series_id: int, path_col: str = \"PNGPath\"):\n _dicom_series_subset_df = dicom_df[dicom_df[\"series_id\"]==series_id] \n return sorted(_dicom_series_subset_df[path_col].tolist(), key=lambda x: int(x.rsplit(\"/\", 1)[-1][:-4]))\n\n\ndef load_dicom_png(png_path: str, n_channels: int = 1, bgr_to_rgb: bool = False):\n \"\"\"\n Load a DICOM image saved as a PNG file.\n\n Args:\n png_path (str): The file path to the PNG image.\n n_channels (int): The number of channels for the output image.\n - 1 for grayscale\n - 3 for RGB color\n\n Returns:\n image (np.ndarray): The loaded image with the specified number of channels.\n\n Raises:\n ValueError: If `n_channels` is not 1 or 3.\n \"\"\"\n # Check if n_channels is valid\n if n_channels not in [1, 3]:\n raise ValueError(\"n_channels must be 1 (grayscale) or 3 (color).\")\n \n # Load the image with the appropriate number of channels\n if n_channels == 1:\n # Load as grayscale\n img = cv2.imread(png_path, cv2.IMREAD_GRAYSCALE)\n elif n_channels == 3:\n # Load as color\n img = cv2.imread(png_path, cv2.IMREAD_COLOR)\n \n # NOTE: \n # - [..., ::-1] to make BGR --> RGB as cv2 loads in BGR as default\n # - Normally this is default/reqd... but in this case it may not matter\n # as the image is likely 1 channel being cast to 3 for annotation purposes\n if bgr_to_rgb:\n img = img[..., ::-1]\n\n # Return the loaded image\n return img\n \n\ndef dicom_df__load_series_images(dicom_df: pd.DataFrame, series_id: int, n_channels: int = 1, bgr_to_rgb: bool = False, path_col: str = \"png_path\", force_channel_dimension: bool = True):\n series_paths = dicom_df__get_series_image_path(dicom_df, series_id, path_col)\n series_img_arr = np.array([load_dicom_png(path, n_channels, bgr_to_rgb) for path in tqdm(series_paths, total=len(series_paths))])\n if len(series_img_arr.shape)==3 and force_channel_dimension:\n series_img_arr = np.expand_dims(series_img_arr, axis=-1)\n return series_img_arr\n\nx = dicom_df__load_series_images(DEMO_DICOM_DF, 3201256954, 1)\nx.shape","metadata":{"execution":{"iopub.status.busy":"2024-06-23T19:01:58.046165Z","iopub.execute_input":"2024-06-23T19:01:58.046546Z","iopub.status.idle":"2024-06-23T19:01:58.793824Z","shell.execute_reply.started":"2024-06-23T19:01:58.046506Z","shell.execute_reply":"2024-06-23T19:01:58.792767Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def create_animation(\n img_stack: np.ndarray, \n study_id: str = \"\",\n series_id: str = \"\",\n imaging_plane: str = \"\",\n figsize: tuple[int, int] = (7, 7), \n interval: int = 83, # 1000//12=83[.333]\n cmap: str = 'bone',\n do_save: bool = False,\n save_dir: str = \"/kaggle/working/animations\"\n) -> animation.FuncAnimation:\n \"\"\"Create an animation from a stack of DICOM images.\n\n Args:\n img_stack (np.ndarray): \n A 3D numpy array representing a stack of DICOM images (shape: [num_frames, height, width, channels]).\n plot_title (str, optional)\n figsize (tuple[int, int], optional): \n A tuple representing the size of the figure (width, height) in inches.\n interval (int, optional): \n Delay between frames in milliseconds. Default is approximately 24 frames per second.\n cmap (str, optional): \n Colormap to use for displaying the images.\n\n Returns:\n animation.FuncAnimation: \n An animation object that can be displayed in a Jupyter notebook or saved to a file.\n \"\"\"\n \n animation_title = f\"STUDY_{study_id}__SERIES_{series_id}__{imaging_plane}_PLANE_VIEW\"\n \n # Create a new figure with the specified size\n fig = plt.figure(figsize=figsize)\n \n # Turn off the axis\n plt.axis('off')\n \n # Display the first image from the stack\n img = plt.imshow(img_stack[0], cmap=cmap)\n plt.title(animation_title.replace(\"__\", \" – \").replace(\"_\", \" \").title(), fontweight=\"bold\")\n def animate_func(i: int) -> list[Any]:\n \"\"\"Update the image for the animation.\n\n Args:\n i (int): The frame index.\n\n Returns:\n list: A list containing the updated image.\n \"\"\"\n img.set_array(img_stack[i])\n return [img]\n plt.close()\n\n # Create and return the animation object\n anim = animation.FuncAnimation(fig, animate_func, frames=len(img_stack), interval=interval)\n \n if do_save:\n if not os.path.isdir(save_dir): \n os.makedirs(save_dir, exist_ok=True)\n anim.save(os.path.join(save_dir, animation_title+\".gif\"), fps=12, writer='imagemagick')\n \n return anim\n\n\n# Create animations for Axial and Sagittal images\ncreate_animation(x, DEMO_STUDY_ID, 3201256954, \"AXIAL\", do_save=True)","metadata":{"execution":{"iopub.status.busy":"2024-06-23T19:01:58.795302Z","iopub.execute_input":"2024-06-23T19:01:58.796120Z","iopub.status.idle":"2024-06-23T19:02:16.308517Z","shell.execute_reply.started":"2024-06-23T19:01:58.796090Z","shell.execute_reply":"2024-06-23T19:02:16.307368Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n\n

Observations

\n\n---\n\n1. Usually 3 view planes\n * Axial T2\n * Sagittal T1\n * Sagittal T2/STIR\n2. In the rare cases where there are fewer than 3 view planes, there will still always be at least 1 Sagittal view and 1 axial view\n3. There are always the same number of images for the pair of Sagittal scans (i.e. the series is complimentary and likely taken at the same time)\n4. Do we need to flip to accomodate...\n* HeadFirstSupine (HFS)\n* FeetFirstSupine (HFS)\n* ImagePatientPosition Z value appears flipped...\n * ```\nPatientPosition\n HFS 118249 \n FFS 28969\n Name: count, dtype: int64\n```","metadata":{}},{"cell_type":"code","source":"for i, (study_id, series_df) in enumerate(train_dicom_df.groupby(\"study_id\")):\n if series_df[\"series_description\"].nunique()!=3:\n print(study_id)\n display(series_df[\"series_description\"].value_counts())\n print()","metadata":{"execution":{"iopub.status.busy":"2024-06-23T19:02:16.309902Z","iopub.execute_input":"2024-06-23T19:02:16.310242Z","iopub.status.idle":"2024-06-23T19:02:16.609137Z","shell.execute_reply.started":"2024-06-23T19:02:16.310213Z","shell.execute_reply":"2024-06-23T19:02:16.607978Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"##### **Class Abstractions**\n\n---\n\nFor a study we know the following:\n- 1 Study is composed of 2-6 Series\n - 3 Studies have 2 Series\n - 1632 Studies have 3 Series \n - 309 Studies have 4 Series\n - 30 Studies have 5 Series\n - 1 Study has 6 Series\n- 1 Study will ALWAYS have at least 1 Axial T2 and 1 Sagittal (T1 or T2) Scan\n - usually though at least 1 count of Axial, 1 count of Sagittal T1 and 1 count of Sagittal T2\n- 1 Study will have 25 labels (probabilities for severity classification)\n - The labels are for multiple conditions (neural foraminal narrowing, subarticular stenosis, and spinal canal stenosis) at five different spinal levels (L1-L2, L2-L3, L3-L4, L4-L5, L5-S1) on both the left and right sides. Each condition and level combination results in a unique label, totaling 25 labels for the study. \n- A Series is either Axial T2, Sagittal T1 or Sagittal T2\n - An Axial T2 scan series is composed of N images/scans:\n - A minimum of 5 images\n - A maximum of 192 images\n - A mean of 34.2 images\n - and 75% of the series will have between 24-44 scans\n - and 99% of the series will have between 14-64 scans\n - A Sagittal T1 scan series is composed of N images/scans **(VERY SIMILAR TO SAGITTAL T2)**:\n - A minimum of 10 images\n - A maximum of 38 images\n - A mean of 17 images\n - and 75% of the series will have between 15-19 scans\n - and 99% of the series will have between 12-25 scans\n - A Sagittal T2 scan series is composed of N images/scans **(VERY SIMILAR TO SAGITTAL T1)**:\n - A minimum of 8 images\n - A maximum of 29 images\n - A mean of 17 images\n - and 75% of the series will have between 15-19 scans\n - and 99% of the series will have between 12-25 scans\n- A single scan (DCM record) has associated metadata, this metadata is useful for providing information about:\n - The image itself (path to the dcm/png and dimensions of the image, pixel spacing, etc.)\n - What the image captures (window size, position of patient/window, slice spacing/thickness, etc.)\n - Rescaling information (RescaleSlope, RescaleType, etc.)\n- As all of the information must be ordered under each relative parent (Instance > Series > Study) and the Study contains the top level labels and predictions, we will make the relevant class structure to allow interaction, manipulation, and visualization of this information accordingly","metadata":{}},{"cell_type":"code","source":"@dataclass\nclass SeriesSlice:\n \"\"\"Represents a single image within a series.\n\n Attributes:\n dcm_path (str): Path to the DICOM file.\n png_path (str): Path to the PNG file.\n study_id (int): ID of the study that contains the series containing this slice.\n series_id (int): ID of the series that contains this slice.\n series_id (int): Description of the series that contains this slice.\n instance_number (int): Instance number of the image.\n img_width (int): Width of the image in pixels.\n img_height (int): Height of the image in pixels.\n pixel_spacing (List[float]): Pixel spacing of the image.\n slice_location (float): Slice location of the image.\n slice_thickness (float): Thickness of the slice.\n slice_spacing (float): Spacing between slices.\n window_center (float): Window center value for the image.\n window_width (float): Window width value for the image.\n RescaleSlope (float): Rescale slope value.\n RescaleType (str): Rescale type.\n SamplesPerPixel (int): Number of samples per pixel.\n primary_axes_with_relative_angle (Any): Primary axes with relative angle.\n \"\"\"\n dcm_path: str\n png_path: str\n study_id: int\n series_id: int\n series_description: str\n instance_number: int\n img_width: int\n img_height: int\n pixel_spacing: list[float]\n slice_location: float\n slice_thickness: float\n slice_spacing: float\n window_center: float\n window_width: float\n RescaleSlope: float\n RescaleType: str\n SamplesPerPixel: int\n primary_axes_with_relative_angle: Any\n\n def load_png(self, n_channels: int = 1, bgr_to_rgb: bool = False) -> np.ndarray:\n \"\"\"Load the image from the PNG file.\n\n Args:\n n_channels (int): Number of channels for the output image.\n bgr_to_rgb (bool): Convert BGR to RGB if True.\n\n Returns:\n np.ndarray: Loaded image.\n \"\"\"\n if n_channels not in [1, 3]:\n raise ValueError(\"n_channels must be 1 (grayscale) or 3 (color).\")\n \n img = cv2.imread(self.png_path, cv2.IMREAD_GRAYSCALE if n_channels == 1 else cv2.IMREAD_COLOR)\n \n if bgr_to_rgb and n_channels == 3:\n img = img[..., ::-1]\n return img\n\n def get_metadata(self) -> dict[str, Any]:\n \"\"\"\n Get metadata associated with the instance.\n\n Returns:\n dict: Metadata dictionary.\n \"\"\"\n return {\n \"dcm_path\": self.dcm_path,\n \"png_path\": self.png_path, \n \"study_id\": self.series_id,\n \"series_id\": self.series_id,\n \"series_description\": self.series_description,\n \"instance_number\": self.instance_number,\n \"img_width\": self.img_width,\n \"img_height\": self.img_height,\n \"pixel_spacing\": self.pixel_spacing,\n \"slice_location\": self.slice_location,\n \"slice_thickness\": self.slice_thickness,\n \"slice_spacing\": self.slice_spacing,\n \"window_center\": self.window_center,\n \"window_width\": self.window_width,\n \"RescaleSlope\": self.RescaleSlope,\n \"RescaleType\": self.RescaleType,\n \"SamplesPerPixel\": self.SamplesPerPixel,\n \"primary_axes_with_relative_angle\": self.primary_axes_with_relative_angle\n }\n\n @classmethod\n def from_row(cls, row: pd.Series) -> 'SeriesSlice':\n \"\"\"\n Create an instance from a DataFrame row.\n\n Args:\n row (pd.Series): DataFrame row.\n\n Returns:\n Instance: Created instance.\n \"\"\"\n return cls(\n dcm_path=row['dcm_path'], \n png_path=row['png_path'], \n study_id=row[\"study_id\"],\n series_id=row[\"series_id\"],\n series_description=row[\"series_description\"],\n instance_number=row['instance_number'], \n img_width=row['img_width'], \n img_height=row['img_height'], \n pixel_spacing=row['pixel_spacing'], \n slice_location=row['slice_location'], \n slice_thickness=row['slice_thickness'], \n slice_spacing=row['slice_spacing'], \n window_center=row['window_center'], \n window_width=row['window_width'], \n RescaleSlope=row['RescaleSlope'], \n RescaleType=row['RescaleType'], \n SamplesPerPixel=row['SamplesPerPixel'], \n primary_axes_with_relative_angle=row['primary_axes_with_relative_angle']\n )\n \n def __repr__(self):\n display(pd.Series(self.get_metadata()))\n return \"\\n\"\n \n\nclass SeriesScan:\n \"\"\"Represents a collection of instances (images/slices) within a study.\n\n Attributes:\n series_df (int): \n The dataframe containing only the subset of the dicom df\n with rows for the series we are examining\n init_series_arr (bool, optional):\n Whether to load the individual images into a numpy array\n\n \"\"\"\n def __init__(self, series_df: pd.DataFrame, init_series_arr: bool = True) -> None:\n self.series_df = series_df[list(SeriesSlice.__annotations__.keys())].sort_values(by=\"instance_number\").reset_index(drop=True)\n self.series_id = self.series_df.series_id.values[0]\n self.series_description = self.series_df.series_description.values[0]\n self.series_slices = [SeriesSlice.from_row(_row) for _, _row in self.series_df.iterrows()]\n self.series_metadata = self.get_metadata()\n self.series_arr = None\n if init_series_arr:\n self.series_arr = self.load_np_arr()\n \n def load_np_arr(\n self, \n bgr_to_rgb: bool = False, \n path_col: str = \"png_path\", \n force_channel_dimension: bool = True\n ) -> np.ndarray:\n \"\"\" TBD \"\"\"\n _img_arr = np.array([_slice.load_png(bgr_to_rgb=bgr_to_rgb) for _slice in self.series_slices])\n if len(_img_arr.shape)==3 and force_channel_dimension:\n _img_arr = np.expand_dims(_img_arr, axis=-1)\n return _img_arr\n \n def create_series_animation(\n self,\n figsize: tuple[int, int] = (7, 7), \n interval: int = 83, # 1000//12=83[.333]\n cmap: str = 'bone',\n do_save: bool = False,\n save_dir: str = \"/kaggle/working/animations\"\n ) -> None:\n \n if self.series_arr is None:\n print(\"\\n... LOADING SERIES NUMPY ARRAY ...\\n\")\n self.series_arr=self.load_np_arr()\n \n display(create_animation(\n self.series_arr, \n study_id=self.series_metadata[\"study_id\"], \n series_id=self.series_id, \n imaging_plane=self.series_description.split(\" \", 1)[0],\n figsize=figsize, interval=interval, cmap=cmap\n ))\n \n def get_metadata(self) -> dict[str, Any]:\n \"\"\"Get metadata associated with the series.\n\n Returns:\n dict: Metadata dictionary.\n \"\"\"\n return {k:v for k,v in self.series_slices[0].get_metadata().items() if k!=\"instance_number\"}\n \n \n @classmethod\n def from_dicom_df(cls, series_id: str | int, dicom_df: pd.DataFrame) -> \"SeriesScan\":\n return cls(dicom_df[dicom_df.series_id==series_id].reset_index(drop=True))\n \n @classmethod\n def from_series_slices(cls, series_slices: list[SeriesSlice]) -> \"SeriesScan\":\n return cls(pd.DataFrame([series_slice.get_metadata() for series_slice in series_slices]))\n \n def __repr__(self):\n display(self.series_df)\n return \"\\n\"\n \n def __getitem__(self, key: int):\n return self.series_slices[key]\n \n def __len__(self):\n return len(self.series_slices)\n\n \nss = SeriesScan.from_dicom_df(702807833, train_dicom_df)\n# ss.create_series_animation()\n# print(len(ss))","metadata":{"execution":{"iopub.status.busy":"2024-06-23T19:02:33.377047Z","iopub.execute_input":"2024-06-23T19:02:33.377453Z","iopub.status.idle":"2024-06-23T19:02:33.521400Z","shell.execute_reply.started":"2024-06-23T19:02:33.377402Z","shell.execute_reply":"2024-06-23T19:02:33.520154Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class PatientStudy:\n \"\"\"Represents a patient's study consisting of multiple series.\n\n Attributes:\n study_df (int): \n The dataframe containing only the subset of the dicom df\n with rows for the study we are examining\n init_study_arr (bool, optional):\n Whether to load the individual images in each series into a numpy array\n force_3_planes (bool, optional):\n Whether to coerce the study so that we always have 3 planes.\n - When more than one series of the same type are found we will\n take the one with more images\n - When less than 1 series of each type/description are found we\n will duplicate the Sagittal scan series and skip this at inference time.\n \"\"\"\n\n def __init__(self, study_df: pd.DataFrame, init_study_arr: bool = True, force_3_planes: bool = True):\n \"\"\"TBD\n \n Args:\n TBD\n \"\"\"\n self.study_df = study_df[list(SeriesSlice.__annotations__.keys())].sort_values(by=[\"series_id\", \"instance_number\"]).reset_index(drop=True)\n self.study_id = self.study_df.study_id.values[0]\n self.series_mapping = self.get_series_mapping([\n SeriesScan.from_dicom_df(series_id, study_df) \n for series_id in study_df.series_id.unique()]\n )\n self.axial_t2_series = self.series_mapping[\"Axial T2\"]\n self.sagittal_t1_series = self.series_mapping[\"Sagittal T1\"]\n self.sagittal_t2_series = self.series_mapping[\"Sagittal T2/STIR\"]\n self.study_metadata = self.get_metadata()\n \n @staticmethod\n def get_series_mapping(series_scans: list[SeriesScan]) -> dict[str, SeriesScan]:\n \"\"\"Process MRI series scans.\n\n Ensures all required types are present, handling missing sagittal scans,\n and selecting the largest scan for each modality.\n\n Args:\n series_scans (list[SeriesScan]): A list of SeriesScan instances to process.\n\n Returns:\n dict[str, SeriesScan]: A dictionary where keys are series descriptions and values\n are the corresponding SeriesScan objects.\n\n Raises:\n ValueError: If the series description is invalid\n ValueError: If the Axial T2 series is missing.\n ValueError: If the both sagittal series are missing.\n\n Example:\n >>> scans = [SeriesScan(\"Axial T2\", [1, 2, 3]), SeriesScan(\"Sagittal T1\", [4, 5])]\n >>> result = process_mri_series(scans)\n >>> print([(k, len(v)) for k, v in result.items()])\n [('Axial T2', 3), ('Sagittal T1', 2), ('Sagittal T2/STIR', 2)]\n \"\"\"\n # Initialize the result dictionary with None values\n result = {series_type: None for series_type in SERIES_TYPES}\n\n # Process each series scan, keeping the largest for each type\n for scan in series_scans:\n if scan.series_description in SERIES_TYPES:\n if result[scan.series_description] is None or len(scan) > len(result[scan.series_description]):\n result[scan.series_description] = scan\n else:\n raise ValueError(f\"\\n...Scan Description ({scan.series_description}) Is Invalid. \"\\\n f\"Should be one of: {SERIES_TYPES}\")\n\n # Check if Axial T2 is present \n if result[\"Axial T2\"] is None:\n raise ValueError(\"\\n... Axial T2 series is missing ...\\n\")\n\n # Check if BOTH Sagittal Scans are missing\n sagittal_scans = [result[\"Sagittal T1\"], result[\"Sagittal T2/STIR\"]]\n if not any(sagittal_scans):\n raise ValueError(\"\\n... Both sagittal scans are missing ...\\n\")\n\n # Handle missing sagittal scans (i.e We have 1 of 2 ... not 2 of 2)\n if not all(sagittal_scans):\n existing_sagittal = next(scan for scan in sagittal_scans if scan is not None)\n for series_type in [\"Sagittal T1\", \"Sagittal T2/STIR\"]:\n if result[series_type] is None:\n result[series_type] = existing_sagittal\n return result\n \n def get_metadata(self) -> dict[str, Any]:\n \"\"\"Get metadata associated with the series.\n\n Returns:\n dict: Metadata dictionary.\n \"\"\"\n return {k:v for k,v in self.axial_t2_series.get_metadata().items() if k not in [\"series_id\", \"series_description\"]}\n\n @classmethod\n def from_dicom_df(cls, study_id: str | int, dicom_df: pd.DataFrame) -> \"SeriesScan\":\n return cls(dicom_df[dicom_df.study_id==study_id].reset_index(drop=True))\n \n def __repr__(self):\n print(\"\\nAXIAL DATA:\\n\")\n display(self.axial_t2_series)\n \n print(\"\\nSAGITTAL T1 DATA:\\n\")\n display(self.sagittal_t1_series)\n \n print(\"\\nSAGITTAL T2 DATA:\\n\")\n display(self.sagittal_t2_series)\n \n return \"\\n\"\n \n def __getitem__(self, key: str):\n if key in self.series_mapping:\n return self.series_mapping[key]\n else:\n raise ValueError(f\"\\n... Key must be one of {SERIES_TYPES} ...\\n\")\n \n ","metadata":{"execution":{"iopub.status.busy":"2024-06-23T19:02:20.460860Z","iopub.execute_input":"2024-06-23T19:02:20.461207Z","iopub.status.idle":"2024-06-23T19:02:20.479261Z","shell.execute_reply.started":"2024-06-23T19:02:20.461178Z","shell.execute_reply":"2024-06-23T19:02:20.478099Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ps = PatientStudy.from_dicom_df(4003253, train_dicom_df)","metadata":{"execution":{"iopub.status.busy":"2024-06-23T19:02:20.480706Z","iopub.execute_input":"2024-06-23T19:02:20.481068Z","iopub.status.idle":"2024-06-23T19:02:21.309566Z","shell.execute_reply.started":"2024-06-23T19:02:20.481039Z","shell.execute_reply":"2024-06-23T19:02:21.308642Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_spine_visualization(axial_df, sagittal_df):\n fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 10))\n fig.suptitle('Spine Visualization: Axial and Sagittal Views', fontsize=16)\n\n # Axial View\n ax1.set_title('Axial View')\n ax1.scatter(axial_df['instance_number'], axial_df['slice_location'], \n c=axial_df['window_center'], cmap='viridis', s=50)\n ax1.set_xlabel('Instance Number')\n ax1.set_ylabel('Slice Location')\n ax1.invert_yaxis() # Invert y-axis to match medical convention\n\n # Add colorbar for axial view\n scatter1 = ax1.scatter(axial_df['instance_number'], axial_df['slice_location'], \n c=axial_df['window_center'], cmap='viridis', s=50)\n plt.colorbar(scatter1, ax=ax1, label='Window Center')\n\n # Sagittal View\n ax2.set_title('Sagittal View')\n ax2.scatter(sagittal_df['instance_number'], sagittal_df['slice_location'], \n c=sagittal_df['window_center'], cmap='plasma', s=50)\n ax2.set_xlabel('Instance Number')\n ax2.set_ylabel('Slice Location')\n\n # Add colorbar for sagittal view\n scatter2 = ax2.scatter(sagittal_df['instance_number'], sagittal_df['slice_location'], \n c=sagittal_df['window_center'], cmap='plasma', s=50)\n plt.colorbar(scatter2, ax=ax2, label='Window Center')\n\n plt.tight_layout()\n plt.show()\n \nplot_spine_visualization(ps.axial_t2_series.series_df, ps.sagittal_t1_series.series_df)","metadata":{"execution":{"iopub.status.busy":"2024-06-23T19:02:21.310799Z","iopub.execute_input":"2024-06-23T19:02:21.311133Z","iopub.status.idle":"2024-06-23T19:02:22.065776Z","shell.execute_reply.started":"2024-06-23T19:02:21.311105Z","shell.execute_reply":"2024-06-23T19:02:22.064709Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_dicom_df[\"slice_location\"]","metadata":{"execution":{"iopub.status.busy":"2024-06-23T19:04:59.413624Z","iopub.execute_input":"2024-06-23T19:04:59.414355Z","iopub.status.idle":"2024-06-23T19:04:59.422633Z","shell.execute_reply.started":"2024-06-23T19:04:59.414320Z","shell.execute_reply":"2024-06-23T19:04:59.421548Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def map_axial_to_sagittal(sagittal_df: pd.DataFrame, axial_df: pd.DataFrame) -> Dict[int, int]:\n \"\"\"Map axial slices to sagittal slice y-coordinates.\n \n Args:\n sagittal_df (pd.DataFrame): DataFrame containing sagittal image data.\n axial_df (pd.DataFrame): DataFrame containing axial image data.\n \n Returns:\n Dict[int, int]: Dictionary mapping sagittal y-coordinates to axial slice indices.\n \"\"\"\n sag_positions = sagittal_df['slice_location'].values\n ax_positions = axial_df['slice_location'].values\n sag_pixel_spacing = sagittal_df['pixel_spacing'].iloc[0][1]\n \n sag_y_axis_to_pixel_space = np.arange(sag_positions.min(), sag_positions.max(), sag_pixel_spacing)\n \n sag_y_coord_to_axial_slice = {}\n for ax_slice_idx, ax_pos in enumerate(ax_positions):\n diffs = np.abs(sag_y_axis_to_pixel_space - ax_pos)\n sag_y_coord = np.argmin(diffs)\n sag_y_coord_to_axial_slice[sag_y_coord] = ax_slice_idx\n \n return sag_y_coord_to_axial_slice\n\ndef map_sagittal_to_axial(sagittal_df: pd.DataFrame, axial_df: pd.DataFrame) -> Dict[int, int]:\n \"\"\"Map sagittal slices to axial slice x-coordinates.\n \n Args:\n sagittal_df (pd.DataFrame): DataFrame containing sagittal image data.\n axial_df (pd.DataFrame): DataFrame containing axial image data.\n \n Returns:\n Dict[int, int]: Dictionary mapping axial x-coordinates to sagittal slice indices.\n \"\"\"\n sag_positions = sagittal_df['slice_location'].values\n ax_pixel_spacing = axial_df['pixel_spacing'].iloc[0][0]\n ax_width = axial_df['img_width'].iloc[0]\n \n ax_x_axis_to_pixel_space = np.arange(0, ax_width * ax_pixel_spacing, ax_pixel_spacing)\n \n ax_x_coord_to_sagittal_slice = {}\n for sag_slice_idx, sag_pos in enumerate(sag_positions):\n diffs = np.abs(ax_x_axis_to_pixel_space - sag_pos)\n ax_x_coord = np.argmin(diffs)\n ax_x_coord_to_sagittal_slice[ax_x_coord] = sag_slice_idx\n \n return ax_x_coord_to_sagittal_slice","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from IPython.display import display\nimport ipywidgets as widgets\nfrom matplotlib.widgets import Slider, RadioButtons\nfrom matplotlib.figure import Figure\nfrom matplotlib.axes import Axes\n\ndef create_enhanced_spine_visualization(\n np_axial: np.ndarray,\n np_sagittal: np.ndarray,\n axial_df: 'pd.DataFrame',\n sagittal_df: 'pd.DataFrame'\n) -> None:\n \"\"\"Create an interactive visualization of spine axial and sagittal views.\n\n This function generates a matplotlib figure with ipywidgets for interactivity\n to explore axial and sagittal views of spine imaging data. It allows\n the user to fix one view and animate the other, with options for slider\n control or automatic animation.\n\n Args:\n np_axial (np.ndarray): 3D numpy array of axial view images.\n np_sagittal (np.ndarray): 3D numpy array of sagittal view images.\n axial_df (pd.DataFrame): DataFrame containing axial metadata.\n sagittal_df (pd.DataFrame): DataFrame containing sagittal metadata.\n \"\"\"\n fig, (ax_sagittal, ax_axial) = plt.subplots(1, 2, figsize=(20, 10))\n fig.suptitle('Enhanced Interactive Spine Visualization', fontsize=16)\n\n # Sagittal View\n ax_sagittal.set_title('Sagittal View')\n sagittal_img = ax_sagittal.imshow(np_sagittal[len(np_sagittal)//2], cmap='gray', aspect='auto')\n ax_sagittal.set_xlabel('Anterior - Posterior')\n ax_sagittal.set_ylabel('Superior - Inferior')\n sagittal_line = ax_sagittal.axhline(y=np_axial.shape[1]//2, color='r', linestyle='--')\n\n # Axial View\n ax_axial.set_title('Axial View')\n axial_img = ax_axial.imshow(np_axial[len(np_axial)//2], cmap='gray', aspect='auto')\n ax_axial.set_xlabel('Left - Right')\n ax_axial.set_ylabel('Anterior - Posterior')\n axial_line = ax_axial.axvline(x=np_sagittal.shape[2]//2, color='r', linestyle='--')\n\n # Add text annotations for patient information\n fig.text(0.01, 0.99, f\"Patient ID: {axial_df['study_id'].iloc[0]}\", ha='left', va='top')\n fig.text(0.01, 0.97, f\"Axial Series: {axial_df['series_description'].iloc[0]}\", ha='left', va='top')\n fig.text(0.01, 0.95, f\"Sagittal Series: {sagittal_df['series_description'].iloc[0]}\", ha='left', va='top')\n\n plt.tight_layout()\n\n # Create ipywidgets\n slice_slider = widgets.IntSlider(\n value=0,\n min=0,\n max=max(len(np_axial), len(np_sagittal)) - 1,\n step=1,\n description='Slice:',\n continuous_update=False\n )\n\n radio_buttons = widgets.RadioButtons(\n options=['Animate Axial', 'Animate Sagittal'],\n description='View to Animate:',\n disabled=False\n )\n\n animation_checkbox = widgets.Checkbox(\n value=False,\n description='Animation Mode',\n disabled=False\n )\n\n play_button = widgets.Play(\n value=0,\n min=0,\n max=max(len(np_axial), len(np_sagittal)) - 1,\n step=1,\n interval=200,\n description=\"Press play\",\n disabled=False\n )\n\n widgets.jslink((play_button, 'value'), (slice_slider, 'value'))\n\n def update(change):\n slice_num = slice_slider.value\n \n if radio_buttons.value == 'Animate Axial':\n axial_img.set_array(np_axial[slice_num % len(np_axial)])\n axial_line.set_xdata(np_sagittal.shape[2]//2)\n sagittal_line.set_ydata(slice_num % np_axial.shape[0])\n ax_axial.set_title(f'Axial View (Slice {slice_num + 1}/{len(np_axial)})')\n ax_sagittal.set_title('Sagittal View (Fixed)')\n else:\n sagittal_img.set_array(np_sagittal[slice_num % len(np_sagittal)])\n sagittal_line.set_ydata(np_axial.shape[1]//2)\n axial_line.set_xdata(slice_num % np_sagittal.shape[0])\n ax_sagittal.set_title(f'Sagittal View (Slice {slice_num + 1}/{len(np_sagittal)})')\n ax_axial.set_title('Axial View (Fixed)')\n \n fig.canvas.draw_idle()\n\n def on_radio_change(change):\n if change['new'] == 'Animate Axial':\n sagittal_img.set_array(np_sagittal[len(np_sagittal)//2])\n slice_slider.max = len(np_axial) - 1\n play_button.max = len(np_axial) - 1\n else:\n axial_img.set_array(np_axial[len(np_axial)//2])\n slice_slider.max = len(np_sagittal) - 1\n play_button.max = len(np_sagittal) - 1\n slice_slider.value = 0\n play_button.value = 0\n update({'new': slice_slider.value})\n\n def on_animation_mode_change(change):\n slice_slider.disabled = change['new']\n play_button.disabled = not change['new']\n\n slice_slider.observe(update, names='value')\n radio_buttons.observe(on_radio_change, names='value')\n animation_checkbox.observe(on_animation_mode_change, names='value')\n\n # Display the widgets and the plot\n display(widgets.VBox([radio_buttons, animation_checkbox, \n widgets.HBox([slice_slider, play_button])]))\n plt.show()\n\ncreate_enhanced_spine_visualization(ps.axial_t2_series.series_arr, ps.sagittal_t1_series.series_arr, ps.axial_t2_series.series_df, ps.sagittal_t1_series.series_df)","metadata":{"execution":{"iopub.status.busy":"2024-06-23T19:04:11.678831Z","iopub.execute_input":"2024-06-23T19:04:11.679294Z","iopub.status.idle":"2024-06-23T19:04:13.147201Z","shell.execute_reply.started":"2024-06-23T19:04:11.679251Z","shell.execute_reply":"2024-06-23T19:04:13.146018Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"_TYPES","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"execution":{"iopub.status.busy":"2024-06-23T19:16:06.972240Z","iopub.execute_input":"2024-06-23T19:16:06.972643Z","iopub.status.idle":"2024-06-23T19:16:06.979170Z","shell.execute_reply.started":"2024-06-23T19:16:06.972612Z","shell.execute_reply":"2024-06-23T19:16:06.978041Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ss.get_metadata()[\"pixel_spacing\"]","metadata":{"execution":{"iopub.status.busy":"2024-06-23T19:17:04.769241Z","iopub.execute_input":"2024-06-23T19:17:04.770025Z","iopub.status.idle":"2024-06-23T19:17:04.776226Z","shell.execute_reply.started":"2024-06-23T19:17:04.769989Z","shell.execute_reply":"2024-06-23T19:17:04.775165Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_series_ipp_arr(series_df: pd.DataFrame, modality: str = \"Axial T2\", reverse_sort: bool = False) -> np.ndarray:\n \"\"\"Extracts and sorts ImagePositionPatient values from a DataFrame.\n\n Args:\n series_df (pd.DataFrame): DataFrame containing the ImagePositionPatient values.\n modality (str): The modality of the series\n reverse_sort (bool): Flag indicating whether to reverse the sort order. Defaults to False.\n\n Returns:\n np.ndarray: Sorted array of ImagePositionPatient values.\n \"\"\"\n \n ipp_xyz_map = {\"sagittal\": 0, \"coronal\": 1, \"axial\": 2}\n ipp_idx = ipp_xyz_map[modality.split(\" \", 1)[0].lower()]\n \n # Extract ImagePositionPatient values and convert them to a NumPy array of floats\n ipp_arr = np.asarray([xyz[ipp_idx] for xyz in series_df[\"ImagePositionPatient\"].tolist()], dtype=np.float32)\n\n # Get the sorting indices based on the positions array and apply to the ipp array\n return ipp_arr[np.argsort(-ipp_arr if reverse_sort else ipp_arr)]\n \n\ndef load_dicom_stack(dicom_folder: str, plane: str, reverse_sort: bool = False) -> dict[str, np.ndarray]:\n \"\"\"Load a stack of DICOM images and convert them to 8-bit.\n \n Args:\n dicom_folder (str): \n Path to the folder containing DICOM files.\n plane (str): \n The plane of the images ('sagittal', 'coronal', 'axial').\n reverse_sort (bool, optional): \n Whether to reverse sort the images.\n \n Returns:\n A dictionary containing the image array, positions, and pixel spacing.\n \"\"\"\n dicom_files = glob(os.path.join(dicom_folder, \"*.dcm\"))\n dicoms = [pydicom.dcmread(f) for f in dicom_files]\n plane_index = {\"sagittal\": 0, \"coronal\": 1, \"axial\": 2}[plane.lower()]\n positions = np.asarray([float(d.ImagePositionPatient[plane_index]) for d in dicoms])\n idx = np.argsort(-positions if reverse_sort else positions)\n ipp = np.asarray([d.ImagePositionPatient for d in dicoms]).astype(\"float\")[idx]\n array = np.stack([d.pixel_array.astype(\"float32\") for d in dicoms])\n array = array[idx]\n return {\n \"array\": convert_to_8bit(array),\n \"positions\": ipp,\n \"pixel_spacing\": np.asarray(dicoms[0].PixelSpacing).astype(\"float\")\n }\n\n\ndef plot_dicom_images(study: pd.DataFrame, image_dir: str) -> None:\n \"\"\"Plot DICOM images for a given study.\n \n Args:\n study (pd.DataFrame): \n DataFrame containing study information.\n image_dir (str): \n Directory containing DICOM images.\n \"\"\"\n for row in study.itertuples():\n if row.series_description == \"Sagittal T2/STIR\":\n sag_t2 = load_dicom_stack(os.path.join(image_dir, str(row.study_id), str(row.series_id)), plane=\"sagittal\")\n elif row.series_description == \"Sagittal T1\":\n sag_t1 = load_dicom_stack(os.path.join(image_dir, str(row.study_id), str(row.series_id)), plane=\"sagittal\")\n elif row.series_description == \"Axial T2\":\n ax_t2 = load_dicom_stack(os.path.join(image_dir, str(row.study_id), str(row.series_id)), plane=\"axial\", reverse_sort=True)\n \n plt.figure(figsize=(14, 5))\n plt.subplot(1, 3, 1)\n plt.imshow(sag_t2[\"array\"][len(sag_t2[\"array\"]) // 2], cmap=\"gray\")\n plt.subplot(1, 3, 2)\n plt.imshow(sag_t1[\"array\"][len(sag_t1[\"array\"]) // 2], cmap=\"gray\")\n plt.subplot(1, 3, 3)\n plt.imshow(ax_t2[\"array\"][len(ax_t2[\"array\"]) // 2], cmap=\"gray\")\n plt.show()\n\n\ndef map_axial_to_sagittal(sag_t2: dict[str, np.ndarray], ax_t2: dict[str, np.ndarray], plot_all_slices: bool = True) -> None:\n \"\"\"Map axial slices to sagittal slices and plot the results.\n \n Args:\n sag_t2 (dict[str, np.ndarray]): \n Dictionary containing sagittal T2 image data.\n ax_t2 (dict[str, np.ndarray]): \n Dictionary containing axial T2 image data.\n plot_all_slices (bool, optional):\n Whether to show each individual side-by-side plot.\n \"\"\"\n top_left_hand_corner_sag_t2 = sag_t2[\"positions\"][len(sag_t2[\"array\"]) // 2]\n sag_y_axis_to_pixel_space = [top_left_hand_corner_sag_t2[2]]\n while len(sag_y_axis_to_pixel_space) < sag_t2[\"array\"].shape[1]:\n sag_y_axis_to_pixel_space.append(sag_y_axis_to_pixel_space[-1] - sag_t2[\"pixel_spacing\"][1])\n \n sag_y_coord_to_axial_slice = {}\n for ax_t2_slice, ax_t2_pos in zip(ax_t2[\"array\"], ax_t2[\"positions\"]):\n diffs = np.abs(np.asarray(sag_y_axis_to_pixel_space) - ax_t2_pos[2])\n sag_y_coord = np.argmin(diffs)\n sag_y_coord_to_axial_slice[sag_y_coord] = ax_t2_slice\n \n sag_midline_slice = sag_t2[\"array\"][len(sag_t2[\"array\"]) // 2]\n plt.imshow(sag_midline_slice, cmap=\"gray\")\n for k in [*sag_y_coord_to_axial_slice]:\n plt.axhline(y=k, color=\"red\", linestyle=\"--\")\n plt.show()\n \n if plot_all_slices:\n for k, v in sag_y_coord_to_axial_slice.items():\n plt.figure(figsize=(10, 4))\n plt.subplot(1, 2, 1)\n plt.imshow(sag_midline_slice, cmap=\"gray\")\n plt.axhline(y=k, color=\"red\", linestyle=\"--\")\n plt.subplot(1, 2, 2)\n plt.imshow(v, cmap=\"gray\")\n plt.show()\n \n\n \n_df = pd.read_csv(\"/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_series_descriptions.csv\")\n_image_dir = \"/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_images/\"\nstudy = _df.loc[_df.study_id == _df.study_id.iloc[10]]\nplot_dicom_images(study, _image_dir)\nax_t2 = sag_t1 = sag_t2 = None \nfor row in study.itertuples():\n if row.series_description == \"Sagittal T2/STIR\":\n sag_t2 = load_dicom_stack(os.path.join(_image_dir, str(row.study_id), str(row.series_id)), plane=\"sagittal\")\n elif row.series_description == \"Sagittal T1\":\n sag_t1 = load_dicom_stack(os.path.join(_image_dir, str(row.study_id), str(row.series_id)), plane=\"sagittal\")\n elif row.series_description == \"Axial T2\":\n ax_t2 = load_dicom_stack(os.path.join(_image_dir, str(row.study_id), str(row.series_id)), plane=\"axial\", reverse_sort=True)\n \n if sag_t2 and ax_t2:\n break\n\n# Assuming the variables sag_t2 and ax_t2 are defined in the study\nmap_axial_to_sagittal(sag_t2, ax_t2)\n\ndel _df, _image_dir, sag_t1, sag_t2, ax_t2, study; gc.collect(); gc.collect();","metadata":{},"execution_count":null,"outputs":[]}]}