{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.12.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## Competition Overview\nThe knee is the most commonly injured and imaged joint in the body, however, the ways in which radiologists interpret MRI imaging scans differ. The goal of the competition is to develop ML models that detect clinically important knee abnormalities, which could provide useful decision support tools to radiologists in practice.\n\n## Competition Training Data & Prior Inventory Findings\n- Basic data inventory and EDA in our notebooks [RSNA_Data_Inventory](https://www.kaggle.com/code/joshuaziel/rsna-data-inventory) and [RSNA_EDA_Preliminary](https://www.kaggle.com/code/joshuaziel/rsna-eda-preliminary).  Many other great notebooks out there on the same.\n- Images and header data themselves are stored in DICOM files associated with a study and various series; all studies have multiple series, and at least one for each anatomical plane.\n- Pixel data in the DICOM files need to be extracted and processed in order to be interpreted by a vision model; headers may contain a number of additional useful labels.\n\n## Image Pre-Processing Background\n- Open weights vision models like DinoV2 expect to recieve images in specific formats.  As an example, for Facebook's DINOv2: \n  - The model expects a patch size of 14, so the dimensions of each image need to be consistently multiples of 14 (eg, 224 x 224 or 512 x 512)\n  - The model expects RGB images, whereas MRI data in the DICOM format is typically greyscale\n  - The model expects 8bit images, but the DICOM format stores image data in 8-bit, 16-bit (12-bit values), or 32-bit integers\n- In addition, we have not deeply explored the distribution of image data in the training DICOM set; the expectation is that there is certainly variation in sizing, but also may be variation in metadata tags that could be important\n\n## Goals for Notebook: \n1) Explore more fulsomely the variation in the actual DICOM data to understand what needs to be hardened for in the pipeline.\n2) Develop a basic image pre-processing pipeline and test-run feature extraction using DINOv2.<>\n\n## Acknowledgements\n- Thanks to `stpete_ishii` for posting a great tutorial notebook on feature visualization: [DINOv2 Feature Grid Visualization](https://www.kaggle.com/code/stpeteishii/dinov2-feature-grid-visualization/notebook)\n","metadata":{}},{"cell_type":"markdown","source":"## 1. Basic Imports and Data Loading\n\n### Key Packages\n1) `torch` and `torchvision` for tensor manipulation\n2) `transformers` for loading DINOv2 for inference\n3) `pydicom` for loading and windowing raw DICOM data\n4) `PIL` (`Image`) and `sklearn.decomposition` (`PCA`) for \n   certain image transforms and feature dimensionality reduction\n   for display\n\n","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np \nimport decimal\nimport torch\nimport matplotlib\nimport matplotlib.pyplot as plt\nimport math\nimport sys\nimport random\nfrom pathlib import Path\nfrom pydicom import dcmread\nfrom pydicom.dataset import FileDataset\nfrom pydicom.pixels import apply_voi_lut\nfrom torchvision.transforms import v2\nfrom transformers import AutoModel\nfrom PIL import Image\nfrom sklearn.decomposition import PCA\nfrom plotnine import *\nfrom tqdm import tqdm\n\nrandom.seed(42)\nDATA_LOC = Path(\"/kaggle/input/competitions/rsna-knee-abnormality-detection/train_series\")\nprint(f\"Data location found: {DATA_LOC.is_dir()}\")\n\ntrain_df = pd.read_csv('/kaggle/input/competitions/rsna-knee-abnormality-detection/train.csv')\nprint(f\"Loaded train.csv: {train_df.shape[0]} rows by {train_df.shape[1]} columns.\")\ntrain_series_df = pd.read_csv('/kaggle/input/competitions/rsna-knee-abnormality-detection/train_series.csv')\nprint(f\"Loaded train_series.csv: {train_series_df.shape[0]} rows by {train_series_df.shape[1]} columns.\")\nlabeled_series_df = pd.merge(\n    train_df.dropna(), \n    train_series_df, \n    how = 'left', \n    on = 'StudyInstanceUID'\n)\nprint(f\"{len(labeled_series_df['StudyInstanceUID'].unique())} fully labeled studies contain {len(labeled_series_df['SeriesInstanceUID'].unique())} series.\")\n\nif Path('/kaggle/input/datasets/joshuaziel/rsna-headers-data/all_headers.csv').is_file():\n   all_headers_df = pd.read_csv('/kaggle/input/datasets/joshuaziel/rsna-headers-data/all_headers.csv', low_memory = False)\n   headers_extracted = True\n   print(f'Headers data found and loaded into all_headers_df: {all_headers_df.shape[0]} rows by {all_headers_df.shape[1]} columns.')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-30T09:50:15.573593Z","iopub.execute_input":"2026-08-30T09:50:15.573932Z","iopub.status.idle":"2026-08-30T09:50:15.897191Z","shell.execute_reply.started":"2026-08-30T09:50:15.57389Z","shell.execute_reply":"2026-08-30T09:50:15.896334Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 2. Series Loader Class (`StudySeries`)\n- Given a base data location (see `DATA_LOC`), as well as a series and study UID\n  loads key DICOM header data, and optionally windowed pixel data by applying\n  the built-in `apply_voi_lut` function from pydicom\n  - **Note** this is dependent upon the DICOM files having the required fields\n    for windowing; unclear whether *all* the data actually do\n- Includes methods to optionally return a series tensor stack or list of images, both\n  sorted on the 'InstanceNumber' DICOM field\n  - **Note** this is dependent upon the DICOM files having the required field; \n    again, unclear whether *all* the data actually do\n- Some basic error handling and logic that might help identify whether things\n  expected to be consistent across a series might vary\n","metadata":{}},{"cell_type":"code","source":"WIDTH = 518 # Divisible into patch size 14\nHEIGHT = 518 # Divisible into patch size 14\nMEAN = (0.485, 0.456, 0.406) # ImageNet Data that DINOv2 was trained on\nSTD = (0.229, 0.224, 0.225) # ImageNet Data that DINOv2 was trained on\n\nclass SeriesPathError(Exception):\n    \"\"\"Path to study series does not exist or is not a directory\"\"\"\n\nclass StudySeries:\n    def __init__(self, study: str, series: str, base: Path = DATA_LOC):\n        series_path = base / study / series\n        \n        if series_path.is_dir():\n            self.path: Path = series_path\n            self.StudyInstanceUID: str = study\n            self.SeriesInstanceUID: str = series\n        else: \n            raise SeriesPathError \n        \n        try: \n            series_header_data = self._parse_series_header_data()\n            series_header_keys = ['series_number', 'acquisition_number', 'manufacturer', 'equipment', 'patient_sex',  \n                                  'laterality', 'img_type', 'photometric_interpretation', 'start_image_position', \n                                  'end_image_position', 'image_orientation', 'num_slices', 'slice_thickness', \n                                  'spacing_btwn_slices', 'rows', 'columns', 'pixel_spacing', 'high_bit', \n                                  'bits_allocated', 'bits_stored', 'window_present']         \n            self.series_header_data: dict = {key: series_header_data[key] for key in series_header_keys}\n            self.dicom_headers:list[FileDataset] | None = series_header_data['dicom_headers']\n            self.pixels_loaded: bool = False\n            self.header_load_error: Exception | None = None\n        except Exception as e:\n            self.series_header_data = None\n            self.dicoms_headers = None\n            self.pixels_loaded = False\n            self.header_load_error: Exception | None = e\n    \n    def get_series_images(self):\n        return self._transform_pixels(\n            [apply_voi_lut(x.get('pixel_array'), x) for x in self._load_study_series(read_pixels = True)],\n            return_image = True\n        )  \n    def get_series_tensor_stack(self):\n         return self._transform_pixels(\n            [apply_voi_lut(x.get('pixel_array'), x) for x in self._load_study_series(read_pixels = True)],\n            return_image = False\n        )\n    def get_series_header_info(self):\n        return {\n            'series_path': self.path,\n            'StudyInstanceUID': self.StudyInstanceUID,\n            'SeriesInstanceUID': self.SeriesInstanceUID,\n            **self.series_header_data\n        }\n\n    def _parse_series_header_data(self):\n        # Header_data contains the sorted list of dicom objects for pydicom\n        header_data = self._load_study_series(read_pixels = False)\n        series_length = len(header_data)\n        \n        # Series level info that is expected to be constant \n        series_number = {x.get('SeriesNumber') for x in header_data}\n        acquisition_number = {x.get('AcquisitionNumber') for x in header_data}\n        manufacturer = {x.get('Manufacturer') for x in header_data}\n        equipment = {x.get('ManufacturerModelName') for x in header_data}\n        patient_sex = {x.get('PatientSex') for x in header_data}\n        laterality = {x.get('Laterality') for x in header_data}\n        image_orientation_patient = {tuple(x.get('ImageOrientationPatient')) for x in header_data}\n        img_type = {tuple(x.get('ImageType')) for x in header_data}\n        photometric_interpretation = {x.get('PhotometricInterpretation') for x in header_data}\n        slice_thickness = {x.get('SliceThickness') for x in header_data}\n        spacing_btwn_slices = {x.get('SpacingBetweenSlices') for x in header_data}\n        rows = {x.get('Rows') for x in header_data}\n        columns = {x.get('Columns') for x in header_data}\n        pixel_spacing = {tuple(x.get('PixelSpacing')) for x in header_data}\n        high_bit = {x.get('HighBit') for x in header_data}\n        bits_allocated = {x.get('BitsAllocated') for x in header_data}\n        bits_stored = {x.get('BitsStored') for x in header_data}\n        window_present = self._check_for_window_data(header_data)\n        return {\n            'series_number': series_number.pop() if len(series_number) == 1 else series_number,\n            'acquisition_number': acquisition_number.pop() if len(acquisition_number) ==1 else acquisition_number,\n            'manufacturer': manufacturer.pop() if len(manufacturer) ==1 else manufacturer,\n            'equipment': equipment.pop() if len(equipment) ==1 else equipment,\n            'patient_sex': patient_sex.pop() if len(patient_sex) == 1 else patient_sex,\n            'laterality': laterality.pop() if len(laterality) == 1 else laterality,\n            'img_type': list(img_type.pop()) if len(img_type) ==1 else img_type,\n            'photometric_interpretation': photometric_interpretation.pop() if len(photometric_interpretation) == 1 else photometric_interpretation,\n            'start_image_position': header_data[0].get('ImagePositionPatient'),\n            'end_image_position': header_data[-1].get('ImagePositionPatient'),\n            'image_orientation': list(image_orientation_patient.pop()) if len(image_orientation_patient) == 1 else image_orientation_patient,\n            'num_slices': series_length,\n            'slice_thickness': slice_thickness.pop() if len(slice_thickness) == 1 else slice_thickness,\n            'spacing_btwn_slices': spacing_btwn_slices.pop() if len(spacing_btwn_slices) ==1 else spacing_btwn_slices,\n            'rows': rows.pop() if len(rows) == 1 else rows,\n            'columns': columns.pop() if len(columns) == 1 else columns,\n            'pixel_spacing': pixel_spacing.pop() if len(pixel_spacing) == 1 else pixel_spacing,\n            'high_bit': high_bit.pop() if len(high_bit) == 1 else high_bit,\n            'bits_allocated': bits_allocated.pop() if len(bits_allocated) == 1 else bits_allocated,\n            'bits_stored': bits_stored.pop() if len(bits_stored) == 1 else bits_stored,\n            'window_present': window_present,\n            'dicom_headers': header_data \n        }\n    \n    def _load_study_series(self, read_pixels: bool = False):\n        no_pixels = not read_pixels\n        dcms = self.path.glob(\"*.dcm\")\n        return sorted(\n            [dcmread(x, stop_before_pixels = no_pixels) for x in dcms],\n            key = lambda x: x.get('InstanceNumber')\n        )\n    def _check_for_window_data(self, header_data: list[dict]) -> bool:\n        window_center = len([x.get('WindowCenter') for x in header_data if x.get('WindowWidth') is not None])\n        window_width = len([x.get('WindowWidth') for x in header_data if x.get('WindowWidth') is not None])\n        return window_center == len(header_data) and window_width == len(header_data)\n\n\n    def _transform_pixels(\n        # Note that this probably doesn't adequately deal with orientation - Fix in future!\n        self,\n        pixels: list[np.ndarray],\n        return_image: bool = False,\n    ) -> list[Image.Image] | torch.Tensor:\n        max_val = float(2 ** self.series_header_data['bits_stored'] - 1)\n        base = v2.Compose([\n            v2.ToImage(),\n            v2.Lambda(lambda t: t / max_val),\n            v2.Resize((HEIGHT, WIDTH), antialias=True),\n            v2.Grayscale(num_output_channels=3),\n        ])\n        if return_image:\n            to_pil = v2.ToPILImage()\n            return [to_pil(base(p)) for p in pixels]\n        norm = v2.Normalize(mean=MEAN, std=STD)\n        return torch.stack([norm(base(p)) for p in pixels], dim=0)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-30T09:50:15.898607Z","iopub.execute_input":"2026-08-30T09:50:15.898988Z","iopub.status.idle":"2026-08-30T09:50:15.922673Z","shell.execute_reply.started":"2026-08-30T09:50:15.898947Z","shell.execute_reply":"2026-08-30T09:50:15.921705Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 3. Assemble Dataset of Key DICOM Header Values\n- **Known issues to address from exploration:**\n  - `Laterality` is probably critical, but is missing from a significant number of\n    records. It can likely be determined from `ImageOrientationPatient` alone or \n    in combination with other values; however, `ImageOrientationPatient` does vary\n    trivially across the dataset - largely based on minor rounding differences\n- **Key Exploratory Analyses:**\n  - Confirm photometric interpretation is `MONOCHROME2` for all samples\n  - Distribution of image and volume sizes, orientations\n  - Slice depth, given overall volume size\n  - Image bit depth distribution","metadata":{}},{"cell_type":"code","source":"class LateralityValueError(Exception):\n    \"\"\"Laterality field contains unexepected values or varies across slices\"\"\"\n\ndef _clean_image_orientation(image_orientation: list | set, dp: int | None = None) -> tuple | None:\n    if isinstance(image_orientation, list):\n        return image_orientation\n    if not isinstance(image_orientation, set):\n        raise TypeError(f\"expected tuple or set, got {type(image_orientation)}\")\n    if not image_orientation:\n        return None\n    if dp is None or dp <= 3: \n        if (dp is not None) and (dp <= 3):\n            tqdm.write(f\"Decimal Places Not Set Too Low ({3}): Starting from the Top\")\n        dp = 0\n        for item in image_orientation:\n            dp = max([dp] + [abs(decimal.Decimal(str(x)).as_tuple().exponent) for x in item])\n\n    while len(image_orientation) > 1 and dp > 3:\n        dp -= 1\n        image_orientation = {tuple(round(x, dp) for x in item) for item in image_orientation}\n\n    return image_orientation.pop() if len(image_orientation) == 1 else None\n\ndef _clean_laterality(laterality: str | set | None) -> str | None:\n        if isinstance(laterality, str) and len(laterality) > 0 and laterality in ['R', 'L', 'U', 'B']:\n            return laterality\n        elif isinstance(laterality, str) and len(laterality) > 0 and laterality in ['RIGHT', 'right', 'Right']:\n            return 'R'\n        elif isinstance(laterality, str) and len(laterality) > 0 and laterality in ['LEFT', 'left', 'Left']:\n            return 'L'\n        elif (isinstance(laterality, str) and len(laterality) == 0) or  (laterality is None):\n            return None\n        else:\n            raise LateralityValueError(f'Value of laterality does not conform:{laterality}')\n\ndef format_header_data(study: str, series:str, base = DATA_LOC) -> dict:\n    s = StudySeries(study, series, base).get_series_header_info()\n    try:\n        start_x, start_y, start_z = s.get('start_image_position')\n    except Exception as e:\n        tqdm.write(f\"Unable to load start position for {study} / {series}: {e}\")\n        start_x, start_y, start_z = None, None, None\n    try:\n         end_x, end_y, end_z = s.get('end_image_position')\n    except Exception as e:\n        tqdm.write(f\"Unable to load end position for {study} / {series}: {e}\")\n        end_x, end_y, end_z = None, None, None\n    try:\n        ax, ay, az, bx, by, bz = _clean_image_orientation(s.get('image_orientation'))\n    except Exception as e:\n        tqdm.write(f\"Unable to load image orientation for {study} / {series}: {e}\")\n        ax, ay, az, bx, by, bz = None, None, None, None, None, None\n    try: \n        laterality = _clean_laterality(s.get('laterality'))\n    except Exception as e:\n        tqdm.write(f\"Unable to load laterality for {study} / {series}: {e}\")\n        laterality = None\n    \n    try: \n        pixel_spacing_r, pixel_spacing_c = s.get('pixel_spacing')\n    except Exception as e:\n        tqdm.write(f\"Unable to load pixel_spacing for {study} / {series}: {e}\")\n        pixel_spacing_r, pixel_spacing_c = None, None\n\n    return {\n        'StudyInstanceUID': s['StudyInstanceUID'],\n        'SeriesInstanceUID': s['SeriesInstanceUID'],\n        'series_number': s.get('series_number'),\n        'acquisition_number': s.get('acquisition_number'),\n        'manufacturer': s.get('manufacturer'),\n        'equipment': s.get('equipment'),\n        'patient_sex': s.get('patient_sex'),\n        'laterality': laterality,\n        'img_type': ','.join(s.get('img_type')),\n        'photometric_interpretation': s.get('photometric_interpretation'),\n        'start_x': start_x,\n        'start_y': start_y,\n        'start_z': start_z,\n        'end_x': end_x,\n        'end_y': end_y,\n        'end_z': end_z,\n        'orientation_Ax': ax,\n        'orientation_Ay': ay,\n        'orientation_Az': az,\n        'orientation_Bx': bx,\n        'orientation_By': by,\n        'orientation_Bz': bz,\n        'num_slices': s.get('num_slices'),\n        'slice_thickness': s.get('slice_thickness'),\n        'spacing_btwn_slices': s.get('spacing_btwn_slices'),\n        'rows': s.get('rows'),\n        'columns': s.get('columns'),\n        'pixel_spacing_r': pixel_spacing_r,\n        'pixel_spacing_c': pixel_spacing_c,\n        'high_bit': s.get('high_bit'),\n        'bits_allocated': s.get('bits_allocated'),\n        'bits_stored': s.get('bits_stored'),\n        'window_present': s.get('window_present')\n    }\n\nif not headers_extracted:\n    all_header_data = []\n    studies = [x.name for x in DATA_LOC.iterdir() if x.is_dir()]\n    for study in tqdm(studies, desc=\"Processing Studies\", mininterval=5.0, dynamic_ncols=True):\n        for series in [x.name for x in (DATA_LOC / study).iterdir() if x.is_dir()]:\n             all_header_data.append(format_header_data(study, series))\n\n    all_headers_df = pd.merge(\n        pd.DataFrame(all_header_data), \n        train_series_df,\n        on = [\"StudyInstanceUID\", \"SeriesInstanceUID\"],\n        how = \"left\"\n    )\n    try: \n        all_headers_df.to_csv('/kaggle/working/all_headers.csv', index = False)\n        print('Headers data saved to file at /kaggle/working/all_headers.csv')\n    except Exception as e:\n        print (f'Error saving headers data: {e}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-30T09:50:15.923829Z","iopub.execute_input":"2026-08-30T09:50:15.924632Z","iopub.status.idle":"2026-08-30T09:50:15.943824Z","shell.execute_reply.started":"2026-08-30T09:50:15.924604Z","shell.execute_reply":"2026-08-30T09:50:15.942978Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Where is there missing data\nmissing_df = (all_headers_df\n    .isna()\n    .sum()\n    .rename_axis('column')\n    .reset_index(name='missing_count')\n    .sort_values('missing_count', ascending = True)\n)\nmissing_df['missing_pct'] = (missing_df['missing_count'] / all_headers_df.shape[0]) * 100\nmissing_df['column'] = pd.Categorical(\n    missing_df['column'],\n    categories=missing_df['column'],\n    ordered=True,\n)\n\n(\n    ggplot(missing_df, aes(x='column', y='missing_pct'))\n    + geom_col(fill = 'blue')\n    + scale_y_continuous(limits=(0, 100), expand=(0, 0))\n    + coord_flip()\n    + labs(\n        title = \"Missing Data\",\n        y = \"Records with Missing Data, %\", \n        x = \"Label\" ,\n    )\n    + theme_minimal()\n    + theme(panel_grid_major_y = element_blank())\n    + theme(plot_title = element_text(weight = 'bold') )\n    + theme(figure_size = (6,6))\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-30T09:50:15.94578Z","iopub.execute_input":"2026-08-30T09:50:15.946146Z","iopub.status.idle":"2026-08-30T09:50:16.391437Z","shell.execute_reply.started":"2026-08-30T09:50:15.94612Z","shell.execute_reply":"2026-08-30T09:50:16.390733Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# All MRI data should be tagged with Photometric Interpretation of 'MONOCHROME2\"\nprint(f\"\"\"\nPhotometric Interpretations in the dataset:\n{'\\n'.join(all_headers_df['photometric_interpretation'].unique())}\n\"\"\")\n\n# Is Windowing data present for `apply_voi_lut` for all Dicoms?\nwindows_absent = (\n    all_headers_df[all_headers_df['window_present'] == False]\n    .apply(\n        lambda x: f\"{x['StudyInstanceUID']} - {x['SeriesInstanceUID']}\",\n        axis = 1\n    )\n    ).tolist()\nif len(windows_absent) == 0:\n    print('All DICOM headers contain windowing information.')\nelse:\n    print(f\"\"\"\nWindowing data missing for {len(windows_absent)} records:\n{'\\n'.join(windows_absent)}\n\"\"\")\n\n# What bit depths are in the data?\nbit_depths = all_headers_df['bits_stored'].value_counts()\nbit_depths = [f\"{x}: {bit_depths[x]} series\" for x in bit_depths.index]\nprint(f\"\"\"\nBit depths (bits stored) present in the data:\n{'\\n'.join(bit_depths)}\n\"\"\")\n\n# Are the MRI images all perfect squares (eg, H == W)?\nnot_square = (\n        all_headers_df[all_headers_df['rows'] != all_headers_df['columns']]\n        .apply(\n            lambda x: f\"{x['StudyInstanceUID']} - {x['SeriesInstanceUID']}: rows = {x['rows']}, columns = {x['columns']}\",\n            axis = 1\n        )\n    ).tolist()\nif len(not_square) == 0:\n    print(\"Row and Column tags in DICOM headers are equal for all MRI images.\")\nelse:\n    print(f\"\"\"\nRow and Column tags differ in {len(not_square)} ({len(not_square) / all_headers_df.shape[0] * 100}%) records.  Examples:\n{('\\n'.join(random.sample(not_square, k = 5)))}\n\"\"\")\n\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-30T09:50:59.903551Z","iopub.execute_input":"2026-08-30T09:50:59.904056Z","iopub.status.idle":"2026-08-30T09:50:59.942992Z","shell.execute_reply.started":"2026-08-30T09:50:59.904021Z","shell.execute_reply":"2026-08-30T09:50:59.942079Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# How Variable is Image Depth and Number of Slices\ng = ggplot(all_headers_df)\n\n(\n    g\n    + geom_boxplot(mapping = aes(x = 'factor(Anatomical_Plane)', y = 'num_slices'), fill = 'blue')\n    + labs(\n        title = \"Number of Slices per Series\",\n        y = \"Number of Slices\", \n        x = \"Frequency, n\" ,\n    )\n    + theme_minimal()\n    + theme(panel_grid_major_x = element_blank())\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-30T09:51:07.861174Z","iopub.execute_input":"2026-08-30T09:51:07.861817Z","iopub.status.idle":"2026-08-30T09:51:08.257713Z","shell.execute_reply.started":"2026-08-30T09:51:07.861784Z","shell.execute_reply":"2026-08-30T09:51:08.257104Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"(\n    g\n    + geom_boxplot(mapping = aes(x = 'factor(Anatomical_Plane)', y = 'slice_thickness'), fill = 'blue')\n    + labs(\n        title = \"Slice Depth by Anatomical Plane\",\n        y = \"Slice Depth, mm\", \n        x = \"Frequency, n\" ,\n    )\n    + theme_minimal()\n    + theme(panel_grid_major_x = element_blank())\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-30T09:51:10.935569Z","iopub.execute_input":"2026-08-30T09:51:10.936058Z","iopub.status.idle":"2026-08-30T09:51:11.393634Z","shell.execute_reply.started":"2026-08-30T09:51:10.936027Z","shell.execute_reply":"2026-08-30T09:51:11.392934Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Conclusions\n- Limited variables missing data from selected DICOM tags\n  - Laterality may be most important and needs to be calculated from ImageOrientationPatient tag data\n- There is a fair amount of pre-processing work to do:\n  - Existing transforms already deal with bit depth; all train series have bit depths of 12 or 16 bits stored\n  - There are a lot of 'non-square' series; depending on the model these may need to be cropped\n  - Image depth varies substantially!","metadata":{}},{"cell_type":"markdown","source":"## 4. Extract and Visualize Features from a Random Series\n- Picking a sagittal series from the labeled data mainly \n  because, not being a radiologist, that plane is most intuitive \n  for me\n- DINOv2 for this initial feature extraction\n  - This means that images should be 8-bit, compatible with a 14 x 14 \n    patch size, RGB (with channels normalized to ImageNet parameters);\n    the methods in the `StudySeries` do this\n- Overlays available for patch means, L2 norms, and PCA (top 3 components)\n  - Primarily focusing on PCA (fit globally on all patches and normalized\n    to global parameters).  \n  - Top 3 parameters pulled for easy visualization\n    on overlays as RGB channels\n  - **Note** as part of a real pipeline, PCA from `sklearn` should be \n    replaced with math done on the tensor with `torch`\n  ","metadata":{}},{"cell_type":"code","source":"# Random Series Selection\n\nrand_study = (\n    labeled_series_df[labeled_series_df['Anatomical_Plane'] == 'Sagittal']\n    .drop(columns = labeled_series_df.columns.difference(['StudyInstanceUID']))\n    .drop_duplicates()\n    .sample(n=1, random_state = 42).iloc[0,0]\n)\nrand_series = (\n    labeled_series_df[(\n        (labeled_series_df['Anatomical_Plane'] == 'Sagittal') &\n        (labeled_series_df['StudyInstanceUID'] == rand_study)\n    )]\n    .drop(columns = labeled_series_df.columns.difference(['SeriesInstanceUID']))\n    .drop_duplicates()\n    .sample(n=1, random_state = 42).iloc[0,0]\n)\n\nrandom_series = StudySeries(rand_study, rand_series)\nprint(f\"\"\"\nSelected random study: {rand_study}\\n\nSelected ransdom series: {rand_series}\\n\nLocation: {random_series.path}\n\"\"\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-30T09:51:15.179124Z","iopub.execute_input":"2026-08-30T09:51:15.179762Z","iopub.status.idle":"2026-08-30T09:51:15.488551Z","shell.execute_reply.started":"2026-08-30T09:51:15.179731Z","shell.execute_reply":"2026-08-30T09:51:15.487767Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# GPU device setting, model Loading, and inference\n# Note that if cuda, conversion to float32 might not be necessary\n\nif torch.backends.mps.is_available():\n    device = torch.device(\"mps\")\nelif torch.cuda.is_available():\n    device = torch.device(\"cuda\")\nelif torch.backends.xla.is_available():\n    device = torch.device(\"xla\")\nelse:\n    device = torch.device(\"cpu\")\nprint(f\"Device set to: {device}\")\n\nmodel = AutoModel.from_pretrained('facebook/dinov2-base').to(device)\nrandom_series_tensor = random_series.get_series_tensor_stack().to(torch.float32).to(device)\n\nwith torch.no_grad():\n    outputs = model(random_series_tensor)\n    last_hidden_states = outputs.last_hidden_state","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-30T09:50:16.406173Z","iopub.status.idle":"2026-08-30T09:50:16.406514Z","shell.execute_reply.started":"2026-08-30T09:50:16.406361Z","shell.execute_reply":"2026-08-30T09:50:16.406385Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"## Process outputs and overlay\n\nclass MisshapenInputError(Exception):\n    \"\"\"Patch size is not a perfect square\"\"\"\n\ndef get_patch_norms(patch_embeddings: torch.Tensor) -> list[np.ndarray]:\n    n_patches = patch_embeddings.size(1)\n    side = round(n_patches ** 0.5)\n    if side * side == n_patches:\n        return [x.detach().cpu().numpy() for x in torch.unbind(torch.linalg.vector_norm(patch_embeddings, dim = 2).reshape(-1, side, side), dim = 0)]\n    else:\n        raise MisshapenInputError\n\ndef get_patch_means(patch_embeddings: torch.Tensor) -> list[np.ndarray]:\n    n_patches = patch_embeddings.size(1)\n    side = round(n_patches ** 0.5)\n    if side * side == n_patches:\n        return [x.detach().cpu().numpy() for x in torch.unbind(patch_embeddings.mean(dim = 2).reshape(-1, side, side), dim = 0)]\n    else:\n        raise MisshapenInputError\n\ndef get_patch_pcas(patch_embeddings: torch.Tensor, n_components: int = 3) -> list[np.ndarray]:\n    n_images, n_patches, n_channels = patch_embeddings.size()\n    side = round(n_patches ** 0.5)\n    \n    if side * side == n_patches:\n        all_patches = patch_embeddings.reshape(-1, n_channels).detach().cpu().numpy()\n        \n        pca = PCA(n_components)\n        pca.fit(all_patches)\n        \n        # 3. Transform each slice using the locked-in master components\n        patch_pcas = []\n        for x in torch.unbind(patch_embeddings, dim=0):\n            patch_pcas.append(pca.transform(x.detach().cpu().numpy()))\n            \n        # 4. Perform series-wide global normalization\n        stacked_pcas = np.stack(patch_pcas)\n        global_min = stacked_pcas.min()\n        global_max = stacked_pcas.max()\n        denom = (global_max - global_min) if (global_max - global_min) != 0 else 1.0\n        \n        return [\n            ((x - global_min) / denom).reshape(side, side, -1) \n            for x in patch_pcas\n        ]\n    else: \n        raise MisshapenInputError\n\ndef get_overlay(study: str, series: str, series_images, embedding_parameter, alpha: float = 0.5) -> None:\n    num_images = len(series_images)\n    ncols = round(num_images ** 0.5)\n    nrows = math.ceil(num_images / ncols)\n    \n    fig, axes = plt.subplots(nrows, ncols, figsize=(ncols * 3, nrows * 3))\n    \n    axes_flat = np.atleast_1d(axes).flat\n\n    for i, (img_base, img_embed) in enumerate(zip(series_images, embedding_parameter)):\n        ax = axes_flat[i]\n\n        ax.imshow(np.asarray(img_base)) \n        \n        ax.imshow(\n            img_embed, \n            alpha=alpha,                  \n            interpolation=\"bicubic\",   \n            extent=[0, 518, 518, 0]     \n        )\n        ax.axis(\"off\") \n        \n    for j in range(num_images, len(axes_flat)):\n        axes_flat[j].axis(\"off\")\n        \n    plt.suptitle(\n        f\"Study: {study}\\nSeries: {series}\", \n        fontsize=16, \n        y=1.01, \n        va=\"bottom\", \n        ha=\"right\"\n    )\n    plt.tight_layout()\n    plt.show()\n\n\n\nptch_embeddings = last_hidden_states[:, 1:, :] # Remove the CLS token\nget_overlay(rand_study, rand_series, random_series.get_series_images(), get_patch_pcas(ptch_embeddings), alpha = 0.6)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-30T09:50:16.407837Z","iopub.status.idle":"2026-08-30T09:50:16.40822Z","shell.execute_reply.started":"2026-08-30T09:50:16.408043Z","shell.execute_reply":"2026-08-30T09:50:16.408069Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 5. WOW - THAT IS PRETTY COOL!\n- Even just first 3 principal components can outline the knee and key structures;\n  no reason that it should be limited to 3 - just good for a visual gut check.","metadata":{}}]}