{"metadata":{"kernelspec":{"display_name":"rsna-knee (3.11.14.final.0)","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.11.14"}},"nbformat_minor":4,"nbformat":4,"cells":[{"id":"0bb2f539","cell_type":"markdown","source":"# What is actually inside the DICOMs\n\nThe pixel data in this competition is ~570 GB across roughly 819,000 slices, and many of our\npreprocessing decisions depend on properties we cannot see in `train.csv`:\nhow many millimetres a pixel covers, how thick a slice is, how big the acquired image is, and\nwhich of the plane × fat-suppression combinations a given study actually contains.\n\nNone of that requires reading pixel data. A DICOM header is a few kilobytes at the front of each\nfile, and `pydicom` can read it without touching the pixels. In this notebook, we read those\nheaders for a sample of studies and plot what we find.\n\nThe data description already covers the file layout and what each CSV column means, so we do not\nrepeat that material here. Instead, we focus on what the headers add.\n\nEverything here is descriptive. Nothing is tuned, and no model is trained.\n","metadata":{}},{"id":"80e8d034","cell_type":"code","source":"import random\nfrom pathlib import Path\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport pydicom\nfrom tqdm import tqdm\n\nROOT = Path(\"/kaggle/input/competitions/rsna-knee-abnormality-detection\")\n\nN_STUDIES = 200          # raise for a fuller picture; 200 finishes in a couple of minutes\nrandom.seed(0)\n\n# Header-only read: a few KB per file instead of a few hundred.\nFIELDS = (\"StudyInstanceUID\", \"SeriesInstanceUID\", \"Rows\", \"Columns\", \"SliceThickness\",\n          \"RepetitionTime\", \"EchoTime\", \"ScanningSequence\", \"InversionTime\")\n\ndef read_header(path):\n    header = pydicom.dcmread(path, stop_before_pixels=True)\n    row = {field: getattr(header, field, None) for field in FIELDS}\n    pixel_spacing = getattr(header, \"PixelSpacing\", None)\n    row[\"spacing\"] = float(pixel_spacing[1]) if pixel_spacing else None\n    return row\n\n# Sorted before sampling so the seed picks the same studies whatever order the disk hands them back.\nstudy_dirs = sorted((ROOT / \"train_series\").iterdir())\nstudy_dirs = random.sample(study_dirs, min(N_STUDIES, len(study_dirs)))\n\ndf = pd.DataFrame([read_header(dcm) for study in tqdm(study_dirs)\n                   for dcm in sorted(study.rglob(\"*.dcm\"))])\nfor column in (\"Rows\", \"Columns\", \"SliceThickness\", \"RepetitionTime\", \"EchoTime\",\n               \"InversionTime\", \"spacing\"):\n    df[column] = pd.to_numeric(df[column], errors=\"coerce\")\n\nprint(f\"{df.StudyInstanceUID.nunique():,} studies · {df.SeriesInstanceUID.nunique():,} series \"\n      f\"· {len(df):,} slices\")\ndf.head(3)","metadata":{},"outputs":[],"execution_count":null},{"id":"f2a90c31","cell_type":"markdown","source":"### Joining the series labels on\n\nThe header table above does not include plane or fat suppression — those labels live in\n`train_series.csv`, which the competition also provides for the test split. We join them onto the\nheader table, then pair the two attributes to define a **slot**: one of the six plane-and-fat-\nsuppression combinations used below. The rest of the cell defines plotting helpers, giving the\ncharts a shared visual style rather than adding analysis.\n","metadata":{}},{"id":"1386490c","cell_type":"code","source":"# The competition CSV labels every series by plane, fat suppression and fluid sensitivity, for\n# the test split as well as train. The slot below uses the fat-suppression axis only -- the data\n# description warns the two protocol flags are correlated but not equivalent in every case.\nseries_csv = pd.read_csv(ROOT / \"train_series.csv\")\ndf = df.merge(series_csv, on=[\"StudyInstanceUID\", \"SeriesInstanceUID\"], how=\"left\")\ndf[\"plane\"] = df[\"Anatomical_Plane\"].str.upper().str[:3]\ndf[\"slot\"] = df[\"plane\"] + np.where(df[\"Fat_Suppression\"].fillna(0).astype(bool), \"_FS\", \"_NOFS\")\n\nBLUE = \"#2a78d6\"\n\ndef frame(ax, xlabel=\"\", ylabel=\"\", grid=\"y\"):\n    for side in (\"top\", \"right\"):\n        ax.spines[side].set_visible(False)\n    ax.set_xlabel(xlabel); ax.set_ylabel(ylabel)\n    if grid:\n        ax.grid(axis=grid, alpha=.3); ax.set_axisbelow(True)\n\ndef barh(counts, xlabel, fmt=\"{:,}\"):\n    \"\"\"Bars biggest-first with the value at the tip. Takes a descending Series.\"\"\"\n    values = counts.to_numpy()\n    positions = np.arange(len(values))[::-1]\n    _, ax = plt.subplots(figsize=(8, 4), dpi=110)\n    ax.barh(positions, values, height=.62, color=BLUE)\n    ax.set_yticks(positions, [str(label) for label in counts.index])\n    for position, value in zip(positions, values):\n        ax.text(value + values.max() * .015, position, fmt.format(value), va=\"center\", fontsize=9)\n    ax.set_xticks([]); frame(ax, xlabel, grid=\"\")\n    ax.spines[\"bottom\"].set_visible(False)\n    plt.show()","metadata":{},"outputs":[],"execution_count":null},{"id":"ce933429","cell_type":"markdown","source":"## In-plane resolution is not one number\n\n`PixelSpacing` is how many millimetres of anatomy one pixel represents — the conversion between image units and the real world. It is **not** constant across this dataset, and the spread between the 5th and 95th percentiles is wide, so two images the same size in pixels can show quite different amounts of anatomy.\n","metadata":{}},{"id":"40b6d1c1","cell_type":"code","source":"spacing = df[\"spacing\"].dropna()\nfig, ax = plt.subplots(figsize=(8, 4), dpi=110)\nax.hist(spacing, bins=np.linspace(0, 1.25, 126), color=BLUE, edgecolor=\"none\")\nframe(ax, \"mm per pixel\", \"slices\")\nplt.show()\nspacing.describe(percentiles=[.05, .25, .5, .75, .95]).round(4)","metadata":{},"outputs":[],"execution_count":null},{"id":"cc4130a6","cell_type":"markdown","source":"## Field of view: how much anatomy is in frame\n\nMultiplying the number of columns by the millimetres each one covers gives the field of view: how much anatomy is in frame. We can interpret this chart alongside the one above — spacing varies partly because the field of view varies, and partly because the pixel grid does.\n","metadata":{}},{"id":"1be070db","cell_type":"code","source":"field_of_view = (df[\"Columns\"] * df[\"spacing\"]).dropna()\nfig, ax = plt.subplots(figsize=(8, 4), dpi=110)\nax.hist(field_of_view, bins=np.linspace(80, 340, 131), color=BLUE, edgecolor=\"none\")\nframe(ax, \"in-plane field of view (mm)\", \"slices\")\nplt.show()\nfield_of_view.describe(percentiles=[.05, .25, .5, .75, .95]).round(1)","metadata":{},"outputs":[],"execution_count":null},{"id":"bed573a7","cell_type":"markdown","source":"## Slice thickness clusters, but the clusters are far apart\n\nEach slice is a slab, and `SliceThickness` is how deep that slab is in millimetres. It is set by the acquisition protocol rather than measured from the image, so it takes a handful of distinct values that span a wide range.\n","metadata":{}},{"id":"75f4d930","cell_type":"code","source":"thickness_counts = df[\"SliceThickness\"].round(2).value_counts().head(10)\nthickness_counts.index = [f\"{mm:g} mm\" for mm in thickness_counts.index]\nbarh(thickness_counts, \"slices\")","metadata":{},"outputs":[],"execution_count":null},{"id":"1cc6a9a2","cell_type":"markdown","source":"## Acquisition matrices are not always square\n\nThe matrix is the pixel grid the scanner acquired, rows by columns. The chart shows the ten most common matrices in the sample; several are not square.\n","metadata":{}},{"id":"152988bc","cell_type":"code","source":"matrix = df[\"Rows\"].astype(\"Int64\").astype(str) + \"x\" + df[\"Columns\"].astype(\"Int64\").astype(str)\nbarh(matrix.value_counts().head(10), \"slices\")","metadata":{},"outputs":[],"execution_count":null},{"id":"cefd7d66","cell_type":"markdown","source":"## Most studies contain only a subset of slots\n\nPairing each of the three planes with fat suppression on or off gives six combinations, which we call **slots**. Each study contains only the slots selected during acquisition: almost no study has all six, and the most common slot appears in roughly an order of magnitude more studies than the rarest.\n","metadata":{}},{"id":"0871357a","cell_type":"code","source":"n_studies = df.StudyInstanceUID.nunique()\nstudies_with_slot = df.dropna(subset=[\"slot\"]).groupby(\"slot\").StudyInstanceUID.nunique()\ncoverage_pct = (100 * studies_with_slot / n_studies).sort_values(ascending=False)\nbarh(coverage_pct, \"% of studies\", \"{:.1f}%\")\ncoverage_pct.round(1)","metadata":{},"outputs":[],"execution_count":null},{"id":"90c184f4","cell_type":"markdown","source":"## Series counts vary by study\n\nThe sample contains no fixed number of series per study; each exam reflects its acquisition protocol.\n","metadata":{}},{"id":"00318439","cell_type":"code","source":"series_per_study = (df.groupby(\"StudyInstanceUID\").SeriesInstanceUID.nunique()\n                    .value_counts().sort_index())\nfig, ax = plt.subplots(figsize=(8, 4), dpi=110)\nax.bar(series_per_study.index, series_per_study.to_numpy(), width=.68, color=BLUE)\nframe(ax, \"series in a study\", \"studies\")\nplt.show()","metadata":{},"outputs":[],"execution_count":null},{"id":"73cb0c2a","cell_type":"markdown","source":"## Contrast weighting across imaging planes\n\nThe CSV gives us the plane and the fat-suppression flag. It does not give us the other axis that\nshapes the appearance of an MR image — its **contrast weighting** — which lives in the header.\n\nTwo header fields help characterise it. `RepetitionTime` (TR) is how long the scanner waits between radio pulses,\nand `EchoTime` (TE) is how long it waits before reading the signal back. Their combination decides\nwhich tissues appear bright:\n\n- **T1** — short TR. Fat bright, fluid dark. Good for anatomy.\n- **T2** — long TR, long TE. Fluid bright, which is what effusion, oedema and fluid-filled tears look like.\n- **PD** — long TR, short TE. Intermediate contrast but a clean signal; the everyday knee sequence.\n\nA third field, `ScanningSequence`, records which *family* of acquisition was used, as a two-letter\ncode: `SE` for spin echo, `GR` for gradient echo, `IR` for inversion recovery. Spin echo is the\ncommon case and the one the three rules above describe. We handle the other two separately:\n\n- **GRE**, from `GR`, is gradient echo. Spin echo fires a second, 180° pulse to bring the signal back\n  into phase; gradient echo instead reverses a magnetic field gradient to do it, and tips the\n  magnetisation by less than 90° to begin with. That is what lets it use a very short TR — so a TR\n  threshold on its own would call it T1.\n- **STIR**, from `IR`, is inversion recovery. An extra 180° pulse flips the magnetisation upside down\n  first, and the image is taken after a delay called `InversionTime` (TI), chosen as the moment fat's\n  signal is passing through zero on its way back. Fat contributes nothing at that instant, so it\n  disappears — suppressed by timing rather than by a fat-selective pulse. Its TR and TE therefore do\n  not describe its appearance either.\n\n`IR` on its own only says *some* inversion recovery: a much longer TI nulls fluid instead of fat and\ngives FLAIR. The TI values printed below are what show these are the fat-nulling kind.\n","metadata":{}},{"id":"a1c4e2b7","cell_type":"markdown","source":"### The values those rules read\n\nThe chart plots one point per series, with TR against TE and the two cut points marked. The\ndistribution also explains the choices of 800 and 60 ms. These are conventional round values from\nradiology teaching rather than physical boundaries, but they are not arbitrary for this sample.\nKnee protocols cluster well to either side of 800 ms, with a nearly empty gap between them, so\nmoving that boundary within the gap relabels almost nothing. The separation around 60 ms is less\ndistinct, making the T2/PD boundary less clear; we can think of T2 and PD as neighbours on a\ncontinuum rather than as clean groups.\n\nBelow the chart, we show the inversion times of the `IR` series and the frequency of each raw\n`ScanningSequence` code.\n","metadata":{}},{"id":"b7d3f019","cell_type":"code","source":"series = df.drop_duplicates(\"SeriesInstanceUID\")\n\nfig, ax = plt.subplots(figsize=(8, 4.4), dpi=110)\nax.scatter(series[\"RepetitionTime\"], series[\"EchoTime\"], s=7, alpha=.3, color=BLUE, edgecolor=\"none\")\nax.axvline(800, color=\"#c0392b\", lw=1)\nax.axhline(60, color=\"#c0392b\", lw=1)\nax.set_xscale(\"log\")\nframe(ax, \"TR (ms, log scale)\", \"TE (ms)\", grid=\"both\")\nplt.show()\nprint(\"inversion time (ms) on the IR series:\")\nprint(series.loc[series[\"ScanningSequence\"].astype(str).str.contains(\"IR\"),\n                 \"InversionTime\"].describe().round(1).to_string())\nseries[\"ScanningSequence\"].astype(str).value_counts()","metadata":{},"outputs":[],"execution_count":null},{"id":"c8e15d42","cell_type":"markdown","source":"### Turning those settings into a label\n\nWe read `ScanningSequence` first: `IR` gives STIR and `GR` gives GRE, because both families require\ndifferent handling. We then classify the remaining spin-echo series using TR and TE: a TR below\n800 ms gives T1; at or above 800 ms, a TE of at least 60 ms gives T2 and a lower TE gives PD.\n\nWe leave a series unlabelled if its header omits TR, TE, or the sequence, so the table below does\nnot include every series in the sample.\n","metadata":{}},{"id":"e1bb1f2a","cell_type":"code","source":"# Contrast weighting from the acquisition parameters. Inversion recovery (STIR) and gradient echo\n# both break a plain TR/TE rule -- GRE runs a short TR by design -- so they are matched first.\nsequence = df[\"ScanningSequence\"].astype(\"string\").fillna(\"\").str.upper()\nrepetition_time, echo_time = df[\"RepetitionTime\"], df[\"EchoTime\"]\ndf[\"weighting\"] = np.select(\n    [sequence.str.contains(\"IR\"),\n     sequence.str.contains(\"GR\"),\n     repetition_time < 800,\n     (repetition_time >= 800) & (echo_time >= 60),\n     (repetition_time >= 800) & (echo_time < 60)],\n    [\"STIR\", \"GRE\", \"T1\", \"T2\", \"PD\"], default=None)\n\ncrosstab = pd.crosstab(df[\"weighting\"], df[\"plane\"])\ncounts = crosstab.to_numpy()\nfig, ax = plt.subplots(figsize=(7, 3.6), dpi=110)\nax.imshow(counts, cmap=\"Blues\", aspect=\"auto\")\nax.set_xticks(range(counts.shape[1]), crosstab.columns)\nax.set_yticks(range(counts.shape[0]), crosstab.index)\nfor (row, col), count in np.ndenumerate(counts):\n    ax.text(col, row, f\"{count:,}\", ha=\"center\", va=\"center\", fontsize=9,\n            color=\"white\" if count > counts.max() * .55 else \"black\")\nax.grid(False); plt.show()\ncrosstab","metadata":{},"outputs":[],"execution_count":null},{"id":"fc80d125","cell_type":"markdown","source":"## Where this leaves us\n\nBy opening the files, we learn five things that are not available in `train.csv`:\n\n1. **Pixel spacing varies substantially**, so a fixed pixel count covers a variable amount of knee.\n2. **Slice thickness varies too**, over a handful of protocol-set values.\n3. **Matrices differ, and some are not square.**\n4. **Slot coverage is uneven** — one plane-and-fat-suppression combination is present in nearly\n   every study, another in a small minority. This comparison requires `train_series.csv` alongside the\n   headers, since the plane and fat-suppression flags come from there.\n5. **Contrast weighting is a third axis**, independent of plane and fat suppression, and only the\n   header carries it.\n\nEverything except point 4 comes from the header alone, at a few kilobytes per file — we decode no\npixel data and never read an image.\n","metadata":{}}]}