{"cells":[{"cell_type":"markdown","source":"# Know Your Knee — EDA, DICOM Mining & Baseline\n\n**RSNA Knee Abnormality Detection** — A comprehensive EDA covering label distributions, radiology report mining, NLP label extraction, DICOM metadata analysis, and a 2.5D CNN baseline that produces a submittable CSV.\n\n**What's inside:**\n1. Data Overview & Structure\n2. Label Deep Dive (12 binary targets, imbalance, co-occurrence)\n   - Labeling Bias Warning\n   - Radiology Report Analysis (sample reports, clinical term frequency)\n   - 2b. NLP Label Extraction — expand from 58 to 4,000+ labeled studies\n3. Series & Study Analysis (planes, sequences)\n4. DICOM Metadata Mining (pixel spacing, slice thickness, manufacturers, field strength)\n5. Image Visualization (windowing, normal vs abnormal, montages)\n6. Slice-Level Analysis\n7. Train vs Test Distribution\n8. Full DICOM Header Reference\n9. Key Takeaways & Modeling Implications\n10. Baseline Model & Submission (2.5D ResNet18, trained on NLP-expanded labels)","metadata":{}},{"cell_type":"code","source":"!pip install matplotlib-venn -q\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport matplotlib.gridspec as gridspec\nimport seaborn as sns\nfrom pathlib import Path\nfrom collections import Counter\nimport pydicom\nimport os\nimport warnings\nwarnings.filterwarnings('ignore')\n\n# Style\nsns.set_theme(style=\"whitegrid\", font_scale=1.1)\nplt.rcParams['figure.dpi'] = 120\nplt.rcParams['figure.facecolor'] = 'white'\n\nCOLORS = sns.color_palette(\"husl\", 12)\n\n# Data directory — mounted via competition data source\nDATA_DIR = Path(\"/kaggle/input/competitions/rsna-knee-abnormality-detection\")\n\nif not DATA_DIR.exists():\n    raise FileNotFoundError(\n        \"Competition data not found! Please add 'rsna-knee-abnormality-detection' \"\n        \"as a data source in notebook settings.\"\n    )\n\nprint(f\"Data directory: {DATA_DIR}\")\nprint(f\"Contents: {sorted([f.name for f in DATA_DIR.iterdir()])}\")\n\nTARGETS = [\n    \"ACL\", \"MCL\", \"Medial Meniscus\", \"Lateral Meniscus\",\n    \"Medial OA\", \"Lateral OA\", \"PF OA\", \"Effusion\",\n    \"Synovitis\", \"Baker's\", \"Contusion\", \"Fracture\"\n]\n\nprint(\"Setup complete!\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n# 1. Data Overview\n\nLet's start by loading all the CSV files and understanding the data structure.","metadata":{}},{"cell_type":"code","source":"# Load all CSV files\ntrain_df = pd.read_csv(DATA_DIR / \"train.csv\")\ntest_df = pd.read_csv(DATA_DIR / \"test.csv\")\ntrain_series_df = pd.read_csv(DATA_DIR / \"train_series.csv\")\ntest_series_df = pd.read_csv(DATA_DIR / \"test_series.csv\")\nsample_sub = pd.read_csv(DATA_DIR / \"sample_submission.csv\")\n\nprint(\"=\" * 60)\nprint(\"DATASET SUMMARY\")\nprint(\"=\" * 60)\nprint(f\"\\ntrain.csv:        {train_df.shape[0]:,} rows × {train_df.shape[1]} cols\")\nprint(f\"test.csv:         {test_df.shape[0]:,} rows × {test_df.shape[1]} cols\")\nprint(f\"train_series.csv: {train_series_df.shape[0]:,} rows × {train_series_df.shape[1]} cols\")\nprint(f\"test_series.csv:  {test_series_df.shape[0]:,} rows × {test_series_df.shape[1]} cols\")\nprint(f\"sample_sub.csv:   {sample_sub.shape[0]:,} rows × {sample_sub.shape[1]} cols\")\n\n# Count DICOM files\ntrain_dcm_dir = DATA_DIR / \"train_series\"\ntest_dcm_dir = DATA_DIR / \"test_series\"\n\ntrain_studies = sorted(train_dcm_dir.iterdir()) if train_dcm_dir.exists() else []\ntest_studies = sorted(test_dcm_dir.iterdir()) if test_dcm_dir.exists() else []\n\nprint(f\"\\nTrain studies (directories): {len(train_studies):,}\")\nprint(f\"Test studies (directories):  {len(test_studies):,}\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Peek at train.csv\nprint(\"train.csv — first 5 rows:\")\nprint(\"-\" * 60)\ndisplay(train_df.head())\nprint(f\"\\nColumns: {list(train_df.columns)}\")\nprint(f\"Missing values:\\n{train_df.isnull().sum()}\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Peek at train_series.csv\nprint(\"train_series.csv — first 5 rows:\")\nprint(\"-\" * 60)\ndisplay(train_series_df.head())\nprint(f\"\\nColumns: {list(train_series_df.columns)}\")\nprint(f\"\\nMissing values:\\n{train_series_df.isnull().sum()}\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Sample submission format\nprint(\"sample_submission.csv:\")\nprint(\"-\" * 60)\ndisplay(sample_sub.head(15))\nprint(f\"\\nUnique studies in submission: {sample_sub.shape[0] // len(TARGETS) if len(TARGETS) > 0 else '?'}\")\nprint(f\"Prediction columns: {[c for c in sample_sub.columns if c != 'StudyInstanceUID']}\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n# 2. Label Deep Dive\n\nThe competition asks us to predict **12 knee abnormalities** from MRI scans. Let's understand the label distribution, class imbalance, and multi-label co-occurrence patterns.\n\n### Clinical context for the 12 targets:\n| Category | Abnormality | What it is |\n|----------|------------|------------|\n| **Ligament** | ACL tear | Anterior cruciate ligament — stabilizes rotation |\n| **Ligament** | MCL tear | Medial collateral ligament — stabilizes inner knee |\n| **Meniscus** | Medial Meniscus tear | Inner cartilage cushion |\n| **Meniscus** | Lateral Meniscus tear | Outer cartilage cushion |\n| **Arthritis** | Medial OA | Osteoarthritis on inner compartment |\n| **Arthritis** | Lateral OA | Osteoarthritis on outer compartment |\n| **Arthritis** | PF OA | Patellofemoral osteoarthritis (kneecap) |\n| **Other** | Effusion | Fluid buildup in joint space |\n| **Other** | Synovitis | Inflammation of synovial membrane |\n| **Other** | Baker's cyst | Fluid-filled cyst behind knee |\n| **Trauma** | Contusion | Bone bruise from impact |\n| **Trauma** | Fracture | Bone break |","metadata":{}},{"cell_type":"code","source":"# Identify target columns (everything except StudyInstanceUID and Report)\ntarget_cols = [c for c in train_df.columns if c not in ('StudyInstanceUID', 'Report')]\nprint(f\"Target columns ({len(target_cols)}): {target_cols}\")\n\n# Only labeled studies (many studies have NaN labels — they have reports but no annotations)\nlabeled_df = train_df.dropna(subset=target_cols, how='all').copy()\nprint(f\"\\nLabeled studies: {len(labeled_df):,} / {len(train_df):,} ({len(labeled_df)/len(train_df)*100:.1f}%)\")\nprint(f\"Unlabeled studies (report-only): {len(train_df) - len(labeled_df):,}\")\n\n# Class distribution (labeled subset only)\nfig, axes = plt.subplots(3, 4, figsize=(18, 12))\nfig.suptitle(f\"Class Distribution for Each Abnormality (n={len(labeled_df)} labeled studies)\", \n             fontsize=16, fontweight='bold', y=1.02)\n\nfor idx, (col, ax) in enumerate(zip(target_cols, axes.ravel())):\n    counts = labeled_df[col].value_counts().sort_index()\n    bars = ax.bar(counts.index.astype(str), counts.values, color=[COLORS[idx], '#ff6b6b'], \n                  edgecolor='white', linewidth=1.5)\n    \n    total = counts.sum()\n    for bar, val in zip(bars, counts.values):\n        ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + total*0.01,\n                f'{val}\\n({val/total*100:.1f}%)', ha='center', va='bottom', fontsize=9)\n    \n    ax.set_title(col, fontweight='bold', fontsize=11)\n    ax.set_xlabel('')\n    ax.set_ylabel('Count')\n    ax.set_xticks([0, 1])\n    ax.set_xticklabels(['Normal (0)', 'Abnormal (1)'])\n\nplt.tight_layout()\nplt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Prevalence rate summary — sorted\nprevalence = labeled_df[target_cols].mean().sort_values(ascending=True)\n\nfig, ax = plt.subplots(figsize=(10, 6))\nbars = ax.barh(range(len(prevalence)), prevalence.values, color=sns.color_palette(\"RdYlGn_r\", len(prevalence)))\nax.set_yticks(range(len(prevalence)))\nax.set_yticklabels(prevalence.index, fontsize=11)\nax.set_xlabel(\"Prevalence (fraction of positive cases)\", fontsize=12)\nax.set_title(\"Abnormality Prevalence — How Imbalanced Is Each Target?\", fontsize=14, fontweight='bold')\n\nfor i, (val, name) in enumerate(zip(prevalence.values, prevalence.index)):\n    ax.text(val + 0.005, i, f'{val:.1%}', va='center', fontsize=10, fontweight='bold')\n\nax.axvline(x=0.5, color='gray', linestyle='--', alpha=0.5, label='Balanced (50%)')\nax.legend()\nplt.tight_layout()\nplt.show()\n\nprint(\"\\nImbalance ratio (negative:positive):\")\nfor col in prevalence.index:\n    pos = labeled_df[col].sum()\n    neg = len(labeled_df) - pos\n    print(f\"  {col:20s}: {neg:,} : {pos:,}  (ratio {neg/max(pos,1):.1f}:1)\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Multi-Label Co-occurrence\nKnee injuries rarely happen in isolation. An ACL tear often comes with bone contusions and effusion. Let's visualize which abnormalities tend to co-occur.","metadata":{}},{"cell_type":"code","source":"# Co-occurrence matrix — P(B=1 | A=1) for each pair\nn_targets = len(target_cols)\ncooccurrence = np.zeros((n_targets, n_targets))\n\nfor i, col_a in enumerate(target_cols):\n    for j, col_b in enumerate(target_cols):\n        mask = labeled_df[col_a] == 1\n        if mask.sum() > 0:\n            cooccurrence[i, j] = labeled_df.loc[mask, col_b].mean()\n\nfig, ax = plt.subplots(figsize=(12, 10))\nmask_diag = np.eye(n_targets, dtype=bool)\nsns.heatmap(cooccurrence, xticklabels=target_cols, yticklabels=target_cols,\n            annot=True, fmt='.2f', cmap='YlOrRd', vmin=0, vmax=1,\n            linewidths=0.5, linecolor='white', ax=ax,\n            mask=mask_diag, cbar_kws={'label': 'P(column | row)'})\n\nax.set_title(\"Multi-Label Co-occurrence: P(column abnormality | row abnormality present)\",\n             fontsize=13, fontweight='bold', pad=15)\nax.set_xlabel(\"Given this is positive...\", fontsize=11)\nax.set_ylabel(\"How often is this also positive?\", fontsize=11)\nplt.xticks(rotation=45, ha='right')\nplt.tight_layout()\nplt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Number of abnormalities per study\nlabeled_df['n_abnormalities'] = labeled_df[target_cols].sum(axis=1)\n\nfig, axes = plt.subplots(1, 2, figsize=(16, 5))\n\n# Distribution of number of abnormalities\ncounts = labeled_df['n_abnormalities'].value_counts().sort_index()\naxes[0].bar(counts.index, counts.values, color=sns.color_palette(\"viridis\", len(counts)),\n            edgecolor='white', linewidth=1)\nfor x, y in zip(counts.index, counts.values):\n    axes[0].text(x, y + 5, f'{y}\\n({y/len(labeled_df)*100:.1f}%)', ha='center', fontsize=9)\naxes[0].set_xlabel(\"Number of abnormalities\", fontsize=12)\naxes[0].set_ylabel(\"Number of studies\", fontsize=12)\naxes[0].set_title(\"How many abnormalities per knee?\", fontsize=13, fontweight='bold')\n\n# Normal vs any abnormality\nhas_any = (labeled_df['n_abnormalities'] > 0).sum()\nnormal = (labeled_df['n_abnormalities'] == 0).sum()\naxes[1].pie([normal, has_any], labels=['Completely Normal', 'At Least 1 Abnormality'],\n            autopct='%1.1f%%', colors=['#2ecc71', '#e74c3c'], startangle=90,\n            textprops={'fontsize': 12}, pctdistance=0.6,\n            explode=[0.05, 0.05])\naxes[1].set_title(\"Normal vs. Abnormal Knees\", fontsize=13, fontweight='bold')\n\nplt.tight_layout()\nplt.show()\n\nprint(f\"\\nMean abnormalities per study: {labeled_df['n_abnormalities'].mean():.2f}\")\nprint(f\"Median: {labeled_df['n_abnormalities'].median():.0f}\")\nprint(f\"Max: {labeled_df['n_abnormalities'].max()}\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Correlation between targets (Phi coefficient = Pearson on binary)\ncorr = labeled_df[target_cols].corr()\n\nfig, ax = plt.subplots(figsize=(11, 9))\nmask_upper = np.triu(np.ones_like(corr, dtype=bool), k=1)\nsns.heatmap(corr, annot=True, fmt='.2f', cmap='RdBu_r', center=0,\n            vmin=-0.3, vmax=0.7, linewidths=0.5, linecolor='white',\n            mask=mask_upper, ax=ax, square=True,\n            cbar_kws={'label': 'Pearson Correlation (Phi Coefficient)'})\nax.set_title(\"Label Correlation Matrix (Phi Coefficient)\", fontsize=14, fontweight='bold', pad=15)\nplt.xticks(rotation=45, ha='right')\nplt.tight_layout()\nplt.show()\n\n# Print strongest correlations\nprint(\"\\nTop 10 strongest positive correlations:\")\ncorr_pairs = []\nfor i in range(len(target_cols)):\n    for j in range(i+1, len(target_cols)):\n        corr_pairs.append((target_cols[i], target_cols[j], corr.iloc[i, j]))\ncorr_pairs.sort(key=lambda x: abs(x[2]), reverse=True)\nfor a, b, r in corr_pairs[:10]:\n    print(f\"  {a:20s} ↔ {b:20s}: r = {r:+.3f}\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### ⚠️ Labeling Bias Warning\n\nThe 58 labeled studies average **4.14 abnormalities per knee** — meaning nearly every labeled study has multiple conditions. This is **not representative** of the full training population (which includes normal knees). The gold labels are heavily biased toward pathological cases, likely because annotators prioritized complex studies.\n\n**Implication:** Prevalence rates computed from gold labels alone (shown above) overestimate how common each condition really is. We need the radiology reports to get a better picture.","metadata":{}},{"cell_type":"markdown","source":"### Radiology Reports — What Do They Look Like?\n\nBefore we try to extract labels from reports, let's peek at what we're working with.","metadata":{}},{"cell_type":"code","source":"# ── Sample reports ────────────────────────────────────────────────────────\nreports = train_df['Report'].dropna()\nprint(f\"Studies with reports: {len(reports):,} / {len(train_df):,} ({len(reports)/len(train_df)*100:.1f}%)\")\nprint(f\"Studies without reports: {train_df['Report'].isna().sum():,}\")\n\n# Show one report from a labeled study and one from an unlabeled study\nlabeled_with_report = train_df[train_df[target_cols].notna().any(axis=1) & train_df['Report'].notna()]\nunlabeled_with_report = train_df[train_df[target_cols].isna().all(axis=1) & train_df['Report'].notna()]\n\nprint(f\"\\n{'='*60}\")\nprint(\"SAMPLE REPORT — Labeled study (has gold labels):\")\nprint(f\"{'='*60}\")\nsample_labeled = labeled_with_report.iloc[0]\nprint(f\"StudyID: {sample_labeled['StudyInstanceUID']}\")\nprint(f\"Labels: {[c for c in target_cols if sample_labeled[c] == 1]}\")\nprint(f\"Report:\\n{sample_labeled['Report'][:500]}\")\n\nprint(f\"\\n{'='*60}\")\nprint(\"SAMPLE REPORT — Unlabeled study (report only, no gold labels):\")\nprint(f\"{'='*60}\")\nsample_unlabeled = unlabeled_with_report.iloc[0]\nprint(f\"StudyID: {sample_unlabeled['StudyInstanceUID']}\")\nprint(f\"Report:\\n{sample_unlabeled['Report'][:500]}\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ── Report statistics + clinical term frequency ──────────────────────────\nimport re as re_mod\n\n# Report length distribution\nreport_lengths = reports.str.len()\n\nfig, axes = plt.subplots(1, 2, figsize=(16, 5))\n\naxes[0].hist(report_lengths, bins=50, color='steelblue', edgecolor='white', alpha=0.8)\naxes[0].set_xlabel(\"Report length (characters)\", fontsize=12)\naxes[0].set_ylabel(\"Count\", fontsize=12)\naxes[0].set_title(\"Report Length Distribution\", fontsize=13, fontweight='bold')\naxes[0].axvline(report_lengths.median(), color='red', linestyle='--', \n                label=f'Median: {report_lengths.median():.0f} chars')\naxes[0].legend()\n\n# Clinical term frequency across all reports\nclinical_terms = {\n    'ACL': r'\\bacl\\b|anterior cruciate',\n    'MCL': r'\\bmcl\\b|medial collateral',\n    'meniscus': r'meniscus|meniscal',\n    'tear': r'\\btear\\b|\\btorn\\b',\n    'effusion': r'effusion',\n    'fracture': r'fracture',\n    'osteoarthritis': r'osteoarthr|arthrosis|\\boa\\b',\n    'contusion': r'contusion|bone bruise',\n    'synovitis': r'synovitis|synovial',\n    'edema': r'\\bedema\\b',\n    'intact': r'\\bintact\\b',\n    'normal': r'\\bnormal\\b',\n    'no ': r'\\bno\\b',\n    \"baker's cyst\": r\"baker|popliteal cyst\",\n    'cartilage': r'cartilage|chondr',\n    'ligament': r'ligament',\n}\n\nreports_lower = reports.str.lower()\nterm_counts = {term: reports_lower.str.contains(pattern, regex=True, na=False).sum() \n               for term, pattern in clinical_terms.items()}\nterm_counts = dict(sorted(term_counts.items(), key=lambda x: x[1], reverse=True))\n\ncolors_bar = sns.color_palette(\"viridis\", len(term_counts))\naxes[1].barh(range(len(term_counts)), list(term_counts.values()), color=colors_bar)\naxes[1].set_yticks(range(len(term_counts)))\naxes[1].set_yticklabels(list(term_counts.keys()), fontsize=10)\naxes[1].set_xlabel(\"Number of reports containing term\", fontsize=12)\naxes[1].set_title(\"Clinical Term Frequency Across All Reports\", fontsize=13, fontweight='bold')\naxes[1].invert_yaxis()\n\nfor i, (term, count) in enumerate(term_counts.items()):\n    axes[1].text(count + len(reports)*0.005, i, f'{count:,} ({count/len(reports)*100:.0f}%)', \n                 va='center', fontsize=9)\n\nplt.tight_layout()\nplt.show()\n\nprint(f\"\\nReport length — min: {report_lengths.min()}, max: {report_lengths.max():,}, \"\n      f\"mean: {report_lengths.mean():.0f}, median: {report_lengths.median():.0f}\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n# 2b. Radiology Report Mining — 75x Label Expansion\n\n**The labeling problem:** Only **58 out of 4,407** training studies (1.3%) have gold-standard binary labels. The remaining 4,349 studies have **radiology reports** but no annotations. Training on just 58 studies is why our baseline AUC is ~0.554 (random chance).\n\n**The fix:** Extract labels from the free-text reports using regex pattern matching. Top Kaggle competitors (Pilkwang Kim, 524 votes) use this approach and report **~82% agreement** with gold labels. It's noisy, but 4,000+ weakly-labeled studies beats 58 perfect ones.\n\n**Our approach:**\n1. For each of the 12 targets, define **positive** patterns (e.g., `acl tear`, `acl rupture`) and **negative** patterns (e.g., `no acl tear`, `acl intact`)\n2. Check negation first — if a negative pattern matches, label = 0\n3. If a positive pattern matches (and no negation), label = 1\n4. If no mention at all, label = **NaN** (not 0 — absence of mention ≠ absence of condition)\n5. Validate against the 58 gold labels to measure extraction quality\n6. Merge: gold labels override NLP where available, NLP fills the rest","metadata":{}},{"cell_type":"code","source":"import re\n\n# ── Pattern definitions for each target ──────────────────────────────────\n# Structure: target_name -> (positive_patterns, negative_patterns)\n# Negative patterns are checked FIRST — if matched, label = 0 even if\n# a positive pattern also matches (handles \"no ACL tear\" correctly).\n\nNLP_PATTERNS = {\n    'ACL': {\n        'pos': [\n            r'acl\\s*(tear|rupture|rupt|torn|injur|disrupt|deficien|incomplet|partial|full[- ]?thickness)',\n            r'anterior\\s+cruciate\\s+(ligament\\s+)?(tear|rupture|torn|injur|disrupt)',\n            r'(torn|ruptured|disrupted|injured)\\s+acl',\n            r'acl\\s+(is\\s+)?(absent|not\\s+intact|not\\s+seen|discontinu)',\n        ],\n        'neg': [\n            r'(no|without|negative\\s+for|absent)\\s+[^.]{0,30}acl\\s*(tear|rupture|injur)',\n            r'acl\\s+(is\\s+)?(intact|normal|unremarkable|preserved|continuous)',\n            r'(intact|normal|stable)\\s+acl',\n        ],\n    },\n    'MCL': {\n        'pos': [\n            r'mcl\\s*(tear|rupture|rupt|torn|injur|sprain|strain|disrupt|partial|full[- ]?thickness)',\n            r'medial\\s+collateral\\s+(ligament\\s+)?(tear|rupture|torn|injur|sprain)',\n            r'(torn|ruptured|disrupted|injured)\\s+mcl',\n        ],\n        'neg': [\n            r'(no|without|negative\\s+for)\\s+[^.]{0,30}mcl\\s*(tear|rupture|injur|sprain)',\n            r'mcl\\s+(is\\s+)?(intact|normal|unremarkable|preserved)',\n            r'(intact|normal|stable)\\s+mcl',\n        ],\n    },\n    'Medial Meniscus': {\n        'pos': [\n            r'medial\\s+meniscus\\s*.{0,20}(tear|torn|rupture|degenerat|macerat|fray|cleavage|flap|bucket|displaced|extru)',\n            r'(tear|torn|rupture).{0,20}medial\\s+meniscus',\n            r'mm\\s+(tear|torn)',  # common abbreviation\n            r'medial\\s+meniscal\\s+(tear|degenerat|fray|cleavage|flap)',\n        ],\n        'neg': [\n            r'(no|without|negative\\s+for)\\s+[^.]{0,30}medial\\s+meniscus\\s*(tear|rupture)',\n            r'(no|without)\\s+[^.]{0,30}medial\\s+meniscal\\s*(tear)',\n            r'medial\\s+meniscus\\s+(is\\s+)?(intact|normal|unremarkable)',\n        ],\n    },\n    'Lateral Meniscus': {\n        'pos': [\n            r'lateral\\s+meniscus\\s*.{0,20}(tear|torn|rupture|degenerat|macerat|fray|cleavage|flap|bucket|displaced|extru)',\n            r'(tear|torn|rupture).{0,20}lateral\\s+meniscus',\n            r'lateral\\s+meniscal\\s+(tear|degenerat|fray|cleavage|flap)',\n        ],\n        'neg': [\n            r'(no|without|negative\\s+for)\\s+[^.]{0,30}lateral\\s+meniscus\\s*(tear|rupture)',\n            r'(no|without)\\s+[^.]{0,30}lateral\\s+meniscal\\s*(tear)',\n            r'lateral\\s+meniscus\\s+(is\\s+)?(intact|normal|unremarkable)',\n        ],\n    },\n    'Medial OA': {\n        'pos': [\n            r'medial\\s+(compartment\\s+)?(osteoarthr|arthrit|arthrosis|OA|degenerat|narrowing|cartilage\\s+loss)',\n            r'(osteoarthr|arthrit|degenerat).{0,30}medial\\s+(compartment|femoral|tibial)',\n            r'medial\\s+(joint\\s+space|compartment)\\s+(narrowing|loss)',\n        ],\n        'neg': [\n            r'(no|without|negative\\s+for)\\s+[^.]{0,30}medial\\s+(compartment\\s+)?(osteoarthr|arthrit|OA|degenerat|narrowing)',\n            r'medial\\s+compartment\\s+(is\\s+)?(normal|preserved|unremarkable)',\n        ],\n    },\n    'Lateral OA': {\n        'pos': [\n            r'lateral\\s+(compartment\\s+)?(osteoarthr|arthrit|arthrosis|OA|degenerat|narrowing|cartilage\\s+loss)',\n            r'(osteoarthr|arthrit|degenerat).{0,30}lateral\\s+(compartment|femoral|tibial)',\n            r'lateral\\s+(joint\\s+space|compartment)\\s+(narrowing|loss)',\n        ],\n        'neg': [\n            r'(no|without|negative\\s+for)\\s+[^.]{0,30}lateral\\s+(compartment\\s+)?(osteoarthr|arthrit|OA|degenerat|narrowing)',\n            r'lateral\\s+compartment\\s+(is\\s+)?(normal|preserved|unremarkable)',\n        ],\n    },\n    'PF OA': {\n        'pos': [\n            r'patellofemoral\\s*(osteoarthr|arthrit|arthrosis|OA|degenerat|narrowing|cartilage\\s+loss|chondr)',\n            r'(PF|patello[- ]?femoral)\\s+(compartment\\s+)?(OA|osteoarthr|degenerat|arthrit|narrowing)',\n            r'patell\\w*\\s+(cartilage\\s+)?(loss|thinning|degenerat|chondr|fissur)',\n        ],\n        'neg': [\n            r'(no|without|negative\\s+for)\\s+[^.]{0,30}patellofemoral\\s*(osteoarthr|arthrit|OA|degenerat)',\n            r'patellofemoral\\s+(compartment\\s+)?(is\\s+)?(normal|preserved|unremarkable)',\n        ],\n    },\n    'Effusion': {\n        'pos': [\n            r'(joint\\s+)?effusion',\n            r'(suprapatellar|joint)\\s+(pouch\\s+)?(fluid|effusion)',\n            r'fluid\\s+(in|within)\\s+(the\\s+)?(joint|knee|suprapatellar)',\n            r'(moderate|large|small|mild|trace|significant)\\s+effusion',\n        ],\n        'neg': [\n            r'no\\s+(joint\\s+)?effusion',\n            r'(no|without|absent|negative\\s+for)\\s+[^.]{0,15}(effusion|joint\\s+fluid)',\n            r'no\\s+(significant\\s+)?(joint\\s+)?fluid',\n        ],\n    },\n    'Synovitis': {\n        'pos': [\n            r'synovitis',\n            r'synovial\\s+(thicken|inflam|hypertroph|proliferat|enhance)',\n        ],\n        'neg': [\n            r'(no|without|absent|negative\\s+for)\\s+[^.]{0,15}synovitis',\n            r'no\\s+[^.]{0,15}synovial\\s+(thicken|inflam)',\n        ],\n    },\n    \"Baker's\": {\n        'pos': [\n            r\"baker'?s?\\s*(cyst|'s)\",\n            r'popliteal\\s+(cyst|fossa\\s+cyst)',\n            r'(cyst|collection)\\s+(in|within|of)\\s+(the\\s+)?popliteal',\n        ],\n        'neg': [\n            r\"(no|without|absent|negative\\s+for)\\s+[^.]{0,15}baker'?s?\\s*cyst\",\n            r'(no|without)\\s+[^.]{0,15}popliteal\\s+cyst',\n        ],\n    },\n    'Contusion': {\n        'pos': [\n            r'(bone\\s+)?contusion',\n            r'bone\\s+(bruise|marrow\\s+edema|marrow\\s+signal)',\n            r'(contusion|bruise)\\s+(of|in|at|within)',\n            r'marrow\\s+edema\\s+(pattern|signal)',\n        ],\n        'neg': [\n            r'(no|without|absent|negative\\s+for)\\s+[^.]{0,15}(contusion|bone\\s+bruise)',\n            r'no\\s+[^.]{0,20}(marrow\\s+edema|bone\\s+marrow\\s+signal)',\n        ],\n    },\n    'Fracture': {\n        'pos': [\n            r'fracture',\n            r'(tibial\\s+plateau|femoral\\s+condyle|patell\\w*|fibula)\\s+fracture',\n            r'(stress|insufficiency|avulsion|osteochondral)\\s+fracture',\n        ],\n        'neg': [\n            r'(no|without|absent|negative\\s+for)\\s+[^.]{0,15}fracture',\n            r'no\\s+(acute\\s+|visible\\s+|definite\\s+)?fracture',\n            r'fracture\\s+(is\\s+)?(excluded|ruled\\s+out|not\\s+seen|absent)',\n        ],\n    },\n}\n\n\ndef extract_labels_from_report(text):\n    \"\"\"Extract 12 binary labels from a radiology report using regex.\n    \n    Returns dict with target names as keys.\n    Values: 1 (positive), 0 (negative/negated), NaN (not mentioned).\n    \"\"\"\n    if pd.isna(text) or not isinstance(text, str) or len(text.strip()) < 10:\n        return {t: np.nan for t in NLP_PATTERNS}\n    \n    text_lower = text.lower()\n    results = {}\n    \n    for target, patterns in NLP_PATTERNS.items():\n        # Check negation first\n        neg_match = any(re.search(p, text_lower) for p in patterns['neg'])\n        pos_match = any(re.search(p, text_lower) for p in patterns['pos'])\n        \n        if neg_match and not pos_match:\n            results[target] = 0\n        elif pos_match and not neg_match:\n            results[target] = 1\n        elif pos_match and neg_match:\n            # Both matched — negation is deliberate, trust it.\n            # e.g. \"No ACL tear\" triggers both pos (acl tear) and neg (no...acl tear).\n            results[target] = 0\n        else:\n            # No mention at all — NaN, not 0\n            results[target] = np.nan\n    \n    return results\n\n\n# Quick sanity check\ntest_report = \"ACL tear with associated bone contusion. No meniscal tear. Moderate joint effusion.\"\ntest_labels = extract_labels_from_report(test_report)\nprint(\"Sanity check on sample report:\")\nprint(f\"  Report: {test_report!r}\")\nprint(f\"  Extracted labels:\")\nfor k, v in test_labels.items():\n    marker = \"+\" if v == 1 else (\"-\" if v == 0 else \"?\")\n    print(f\"    [{marker}] {k}: {v}\")\nprint(f\"\\nPatterns defined for {len(NLP_PATTERNS)} targets\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ── Apply NLP extraction to ALL training studies ─────────────────────────\nnlp_labels = train_df['Report'].apply(extract_labels_from_report)\nnlp_df = pd.DataFrame(nlp_labels.tolist(), index=train_df.index)\n\n# How many studies got at least one label?\nhas_any_nlp = nlp_df.notna().any(axis=1).sum()\nhas_no_report = train_df['Report'].isna().sum()\n\nprint(f\"Total training studies:         {len(train_df):,}\")\nprint(f\"Studies with reports:            {len(train_df) - has_no_report:,}\")\nprint(f\"Studies with >= 1 NLP label:     {has_any_nlp:,}\")\nprint(f\"Studies with no NLP labels:      {len(train_df) - has_any_nlp:,}\")\nprint(f\"Gold-labeled studies:            {len(labeled_df):,}\")\n\n# ── NLP coverage per target ──────────────────────────────────────────────\nprint(\"\\n\" + \"=\" * 65)\nprint(f\"{'Target':20s} {'Mentioned':>10s} {'Positive':>10s} {'Negative':>10s} {'No mention':>10s}\")\nprint(\"-\" * 65)\nfor col in target_cols:\n    mentioned = nlp_df[col].notna().sum()\n    pos = (nlp_df[col] == 1).sum()\n    neg = (nlp_df[col] == 0).sum()\n    nan_count = nlp_df[col].isna().sum()\n    print(f\"{col:20s} {mentioned:>10,} {pos:>10,} {neg:>10,} {nan_count:>10,}\")\n\n# ── Prevalence comparison: Gold vs NLP ───────────────────────────────────\nfig, ax = plt.subplots(figsize=(12, 6))\n\ngold_prev = labeled_df[target_cols].mean().values\nnlp_prev = nlp_df[target_cols].apply(lambda x: x.dropna().mean()).values\n\nx = np.arange(len(target_cols))\nw = 0.35\nbars1 = ax.bar(x - w/2, gold_prev, w, label=f'Gold labels (n={len(labeled_df)})', color='#3498db', edgecolor='white')\nbars2 = ax.bar(x + w/2, nlp_prev, w, label=f'NLP-extracted (n={has_any_nlp})', color='#e74c3c', edgecolor='white')\n\nax.set_xticks(x)\nax.set_xticklabels(target_cols, rotation=45, ha='right')\nax.set_ylabel(\"Prevalence (fraction positive)\")\nax.set_title(\"Label Prevalence — Gold Standard vs NLP Extraction\", fontsize=14, fontweight='bold')\nax.legend()\nax.set_ylim(0, 1)\n\nplt.tight_layout()\nplt.show()\n\nprint(\"\\nNote: Gold labels are biased toward pathological cases (avg 4.14 abnormalities/study).\")\nprint(\"NLP prevalence is likely closer to the true population distribution.\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ── Validate NLP labels against gold standard ────────────────────────────\nfrom sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score\n\n# Get NLP predictions for the 58 gold-labeled studies\ngold_idx = labeled_df.index\nnlp_for_gold = nlp_df.loc[gold_idx]\n\nprint(\"NLP Label Extraction — Validation Against 58 Gold Labels\")\nprint(\"=\" * 80)\nprint(f\"{'Target':20s} {'Acc':>6s} {'Prec':>6s} {'Recall':>6s} {'F1':>6s}  {'TP':>4s} {'FP':>4s} {'FN':>4s} {'TN':>4s}  {'Cov':>5s}\")\nprint(\"-\" * 80)\n\noverall_correct = 0\noverall_total = 0\n\nfor col in target_cols:\n    gold = labeled_df[col].values\n    nlp = nlp_for_gold[col].values\n    \n    # Only evaluate where NLP made a prediction (not NaN)\n    mask = ~np.isnan(nlp)\n    coverage = mask.sum() / len(gold)\n    \n    if mask.sum() == 0:\n        print(f\"{col:20s} {'—':>6s} {'—':>6s} {'—':>6s} {'—':>6s}  {'—':>4s} {'—':>4s} {'—':>4s} {'—':>4s}  {coverage:>5.0%}\")\n        continue\n    \n    g, n = gold[mask], nlp[mask].astype(int)\n    \n    tp = ((g == 1) & (n == 1)).sum()\n    fp = ((g == 0) & (n == 1)).sum()\n    fn = ((g == 1) & (n == 0)).sum()\n    tn = ((g == 0) & (n == 0)).sum()\n    \n    acc = accuracy_score(g, n)\n    prec = precision_score(g, n, zero_division=0)\n    rec = recall_score(g, n, zero_division=0)\n    f1 = f1_score(g, n, zero_division=0)\n    \n    overall_correct += (tp + tn)\n    overall_total += len(g)\n    \n    print(f\"{col:20s} {acc:>6.1%} {prec:>6.1%} {rec:>6.1%} {f1:>6.1%}  {tp:>4d} {fp:>4d} {fn:>4d} {tn:>4d}  {coverage:>5.0%}\")\n\noverall_acc = overall_correct / overall_total if overall_total > 0 else 0\nprint(\"-\" * 80)\nprint(f\"{'OVERALL':20s} {overall_acc:>6.1%}\")\nprint(f\"\\nCov = coverage (fraction of gold studies where NLP made a prediction)\")\nprint(f\"Low coverage targets may need more pattern variants or multilingual support.\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ── Merge: Gold labels override NLP, NLP fills the rest ──────────────────\n# Start with NLP labels for all studies\nlabeled_df_expanded = train_df[['StudyInstanceUID']].copy()\nfor col in target_cols:\n    labeled_df_expanded[col] = nlp_df[col]\n\n# Override with gold labels where available (gold is ground truth)\nfor col in target_cols:\n    gold_mask = labeled_df.index\n    labeled_df_expanded.loc[gold_mask, col] = labeled_df[col].values\n\n# Drop rows where ALL targets are NaN (no info at all)\nlabeled_df_expanded = labeled_df_expanded.dropna(subset=target_cols, how='all').copy()\n\n# For training, fill remaining NaN with 0 (conservative: assume normal if not mentioned)\n# Keep a record of which were imputed for potential label-smoothing later\nnlp_imputed_mask = labeled_df_expanded[target_cols].isna()\nn_imputed = nlp_imputed_mask.sum().sum()\nlabeled_df_expanded[target_cols] = labeled_df_expanded[target_cols].fillna(0)\n\nprint(\"=\" * 60)\nprint(\"EXPANDED DATASET SUMMARY\")\nprint(\"=\" * 60)\nprint(f\"  Gold-labeled studies:      {len(labeled_df):>6,}\")\nprint(f\"  NLP-labeled studies:       {len(labeled_df_expanded) - len(labeled_df):>6,}\")\nprint(f\"  Total expanded dataset:    {len(labeled_df_expanded):>6,}\")\nprint(f\"  Expansion factor:          {len(labeled_df_expanded)/len(labeled_df):>6.1f}x\")\nprint(f\"  NaN cells filled with 0:   {n_imputed:>6,} ({n_imputed/(len(labeled_df_expanded)*len(target_cols))*100:.1f}%)\")\n\n# Show new prevalence\nprint(f\"\\nExpanded dataset prevalence:\")\nfor col in target_cols:\n    prev = labeled_df_expanded[col].mean()\n    gold_prev_val = labeled_df[col].mean()\n    print(f\"  {col:20s}: {prev:>6.1%}  (was {gold_prev_val:>6.1%} in gold)\")\n\nprint(f\"\\n>>> labeled_df_expanded is ready for training ({len(labeled_df_expanded):,} studies, {len(target_cols)} targets)\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n# 3. Series & Study Analysis\n\nEach patient study contains multiple MRI **series** — different imaging sequences acquired in different anatomical planes. Understanding this structure is key to building the right model architecture.","metadata":{}},{"cell_type":"code","source":"# Series per study\nseries_per_study = train_series_df.groupby('StudyInstanceUID').size()\n\nfig, axes = plt.subplots(1, 2, figsize=(16, 5))\n\n# Histogram of series per study\naxes[0].hist(series_per_study.values, bins=range(1, series_per_study.max()+2),\n             color='steelblue', edgecolor='white', linewidth=1, align='left')\naxes[0].set_xlabel(\"Number of series per study\", fontsize=12)\naxes[0].set_ylabel(\"Number of studies\", fontsize=12)\naxes[0].set_title(\"How many MRI series per patient?\", fontsize=13, fontweight='bold')\naxes[0].axvline(series_per_study.median(), color='red', linestyle='--', label=f'Median: {series_per_study.median():.0f}')\naxes[0].legend()\n\n# Series descriptions — what types of MRI sequences are used?\ndesc_counts = train_series_df['Anatomical_Plane'].value_counts().head(20)\naxes[1].barh(range(len(desc_counts)), desc_counts.values, color='teal')\naxes[1].set_yticks(range(len(desc_counts)))\naxes[1].set_yticklabels(desc_counts.index, fontsize=9)\naxes[1].set_xlabel(\"Count\", fontsize=12)\naxes[1].set_title(\"Anatomical Plane Distribution\", fontsize=13, fontweight='bold')\naxes[1].invert_yaxis()\n\nplt.tight_layout()\nplt.show()\n\nprint(f\"\\nSeries per study — min: {series_per_study.min()}, max: {series_per_study.max()}, \"\n      f\"mean: {series_per_study.mean():.1f}, median: {series_per_study.median():.0f}\")\nprint(f\"Unique series descriptions: {train_series_df['Anatomical_Plane'].nunique()}\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Anatomical plane analysis\n# The Anatomical_Plane column from the CSV already classifies planes\n# For DICOM metadata, we also define a helper to classify from SeriesDescription\ndef classify_plane(desc):\n    \"\"\"Map series/plane description to standard anatomical plane name.\"\"\"\n    if pd.isna(desc):\n        return 'Unknown'\n    desc = str(desc).lower().strip()\n    if desc in ('sagittal', 'sag'):\n        return 'Sagittal'\n    elif desc in ('coronal', 'cor'):\n        return 'Coronal'\n    elif desc in ('axial', 'ax', 'transverse', 'tra'):\n        return 'Axial'\n    elif 'sag' in desc:\n        return 'Sagittal'\n    elif 'cor' in desc:\n        return 'Coronal'\n    elif 'ax' in desc or 'tra' in desc:\n        return 'Axial'\n    else:\n        return desc.title()  # Use as-is for unknown planes\n\ntrain_series_df['plane'] = train_series_df['Anatomical_Plane'].apply(classify_plane)\n\nplane_counts = train_series_df['plane'].value_counts()\nprint(f\"Plane values: {dict(plane_counts)}\")\n\nfig, axes = plt.subplots(1, 2, figsize=(14, 5))\n\n# Plane distribution\ncolors_plane = {'Sagittal': '#3498db', 'Coronal': '#e74c3c', 'Axial': '#2ecc71', 'Other': '#95a5a6', 'Unknown': '#bdc3c7'}\naxes[0].pie(plane_counts.values, labels=plane_counts.index, autopct='%1.1f%%',\n            colors=[colors_plane.get(p, '#95a5a6') for p in plane_counts.index],\n            startangle=90, textprops={'fontsize': 11})\naxes[0].set_title(\"Distribution of Anatomical Planes\", fontsize=13, fontweight='bold')\n\n# Planes per study — does every study have all 3 planes?\nplanes_per_study = train_series_df.groupby('StudyInstanceUID')['plane'].apply(lambda x: set(x))\nplane_combos = planes_per_study.apply(lambda x: ', '.join(sorted(x))).value_counts().head(10)\naxes[1].barh(range(len(plane_combos)), plane_combos.values, color='mediumpurple')\naxes[1].set_yticks(range(len(plane_combos)))\naxes[1].set_yticklabels(plane_combos.index, fontsize=9)\naxes[1].set_xlabel(\"Number of studies\", fontsize=12)\naxes[1].set_title(\"Most common plane combinations per study\", fontsize=13, fontweight='bold')\naxes[1].invert_yaxis()\n\nplt.tight_layout()\nplt.show()\n\n# Studies with all 3 main planes\nhas_all_3 = planes_per_study.apply(lambda x: {'Sagittal', 'Coronal', 'Axial'}.issubset(x)).sum()\nprint(f\"\\nStudies with all 3 planes (Sag+Cor+Ax): {has_all_3:,} / {len(planes_per_study):,} ({has_all_3/len(planes_per_study)*100:.1f}%)\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n# 4. DICOM Metadata Mining\n\nThis is where it gets interesting. DICOM files contain rich metadata beyond just the pixel data. Let's mine this to understand the imaging parameters, which are crucial for preprocessing decisions.\n\nWe'll sample DICOM files to extract:\n- **Pixel Spacing** — physical size of each pixel (mm)\n- **Slice Thickness** — distance between slices (mm)  \n- **Rows/Columns** — image dimensions\n- **Manufacturer** — scanner brand\n- **Magnetic Field Strength** — 1.5T vs 3T\n- **Sequence Name** — pulse sequence type\n- **Window Center/Width** — display windowing","metadata":{}},{"cell_type":"code","source":"%%time\n# Sample DICOM metadata from training set\n# We sample 1 DICOM per series for speed (metadata is consistent within a series)\n\nimport random\nrandom.seed(42)\n\ndcm_metadata = []\nsample_studies = random.sample(train_studies, min(300, len(train_studies)))\n\nfor study_dir in sample_studies:\n    study_id = study_dir.name\n    if not study_dir.is_dir():\n        continue\n    for series_dir in study_dir.iterdir():\n        if not series_dir.is_dir():\n            continue\n        series_id = series_dir.name\n        dcm_files = list(series_dir.glob(\"*.dcm\"))\n        if not dcm_files:\n            continue\n        \n        # Read one DICOM file per series for metadata\n        dcm = pydicom.dcmread(str(dcm_files[0]), stop_before_pixels=True)\n        \n        row = {\n            'StudyInstanceUID': study_id,\n            'SeriesInstanceUID': series_id,\n            'n_slices': len(dcm_files),\n            'rows': getattr(dcm, 'Rows', None),\n            'columns': getattr(dcm, 'Columns', None),\n            'pixel_spacing_r': float(dcm.PixelSpacing[0]) if hasattr(dcm, 'PixelSpacing') else None,\n            'pixel_spacing_c': float(dcm.PixelSpacing[1]) if hasattr(dcm, 'PixelSpacing') else None,\n            'slice_thickness': float(dcm.SliceThickness) if hasattr(dcm, 'SliceThickness') else None,\n            'manufacturer': getattr(dcm, 'Manufacturer', 'Unknown'),\n            'model': getattr(dcm, 'ManufacturerModelName', 'Unknown'),\n            'field_strength': float(dcm.MagneticFieldStrength) if hasattr(dcm, 'MagneticFieldStrength') else None,\n            'sequence_name': getattr(dcm, 'SequenceName', 'Unknown'),\n            'series_description': getattr(dcm, 'SeriesDescription', 'Unknown'),\n            'laterality': getattr(dcm, 'Laterality', getattr(dcm, 'ImageLaterality', 'Unknown')),\n            'bits_stored': getattr(dcm, 'BitsStored', None),\n            'window_center': dcm.WindowCenter if hasattr(dcm, 'WindowCenter') else None,\n            'window_width': dcm.WindowWidth if hasattr(dcm, 'WindowWidth') else None,\n            'photometric': getattr(dcm, 'PhotometricInterpretation', 'Unknown'),\n        }\n        dcm_metadata.append(row)\n\ndcm_df = pd.DataFrame(dcm_metadata)\nprint(f\"Sampled DICOM metadata from {dcm_df['StudyInstanceUID'].nunique()} studies, {len(dcm_df)} series\")\ndisplay(dcm_df.head())","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Image dimensions\nfig, axes = plt.subplots(2, 2, figsize=(16, 12))\n\n# Rows x Columns scatter\nax = axes[0, 0]\nax.scatter(dcm_df['columns'], dcm_df['rows'], alpha=0.3, s=10, c='steelblue')\nax.set_xlabel(\"Columns (width)\", fontsize=12)\nax.set_ylabel(\"Rows (height)\", fontsize=12)\nax.set_title(\"Image Dimensions (Rows × Columns)\", fontsize=13, fontweight='bold')\n\n# Resolution distribution\ndim_counts = dcm_df.apply(lambda r: f\"{int(r['rows'])}×{int(r['columns'])}\" if pd.notna(r['rows']) else 'Unknown', axis=1)\ntop_dims = dim_counts.value_counts().head(10)\nax = axes[0, 1]\nax.barh(range(len(top_dims)), top_dims.values, color='coral')\nax.set_yticks(range(len(top_dims)))\nax.set_yticklabels(top_dims.index)\nax.set_xlabel(\"Count\", fontsize=12)\nax.set_title(\"Top 10 Image Resolutions\", fontsize=13, fontweight='bold')\nax.invert_yaxis()\n\n# Pixel spacing\nax = axes[1, 0]\nvalid_ps = dcm_df.dropna(subset=['pixel_spacing_r'])\nax.hist(valid_ps['pixel_spacing_r'], bins=50, color='seagreen', edgecolor='white', alpha=0.8)\nax.set_xlabel(\"Pixel Spacing (mm)\", fontsize=12)\nax.set_ylabel(\"Count\", fontsize=12)\nax.set_title(\"Pixel Spacing Distribution\", fontsize=13, fontweight='bold')\nax.axvline(valid_ps['pixel_spacing_r'].median(), color='red', linestyle='--', \n           label=f\"Median: {valid_ps['pixel_spacing_r'].median():.3f} mm\")\nax.legend()\n\n# Slice thickness\nax = axes[1, 1]\nvalid_st = dcm_df.dropna(subset=['slice_thickness'])\nax.hist(valid_st['slice_thickness'], bins=50, color='mediumpurple', edgecolor='white', alpha=0.8)\nax.set_xlabel(\"Slice Thickness (mm)\", fontsize=12)\nax.set_ylabel(\"Count\", fontsize=12)\nax.set_title(\"Slice Thickness Distribution\", fontsize=13, fontweight='bold')\nax.axvline(valid_st['slice_thickness'].median(), color='red', linestyle='--',\n           label=f\"Median: {valid_st['slice_thickness'].median():.1f} mm\")\nax.legend()\n\nplt.tight_layout()\nplt.show()\n\nprint(f\"\\nPixel spacing range: {valid_ps['pixel_spacing_r'].min():.3f} – {valid_ps['pixel_spacing_r'].max():.3f} mm\")\nprint(f\"Slice thickness range: {valid_st['slice_thickness'].min():.1f} – {valid_st['slice_thickness'].max():.1f} mm\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Manufacturer and scanner analysis\nfig, axes = plt.subplots(1, 3, figsize=(20, 5))\n\n# Manufacturers\nmfr_counts = dcm_df['manufacturer'].value_counts().head(10)\naxes[0].barh(range(len(mfr_counts)), mfr_counts.values, color='steelblue')\naxes[0].set_yticks(range(len(mfr_counts)))\naxes[0].set_yticklabels(mfr_counts.index, fontsize=9)\naxes[0].set_xlabel(\"Count\")\naxes[0].set_title(\"Scanner Manufacturers\", fontsize=13, fontweight='bold')\naxes[0].invert_yaxis()\n\n# Field strength\nfs = dcm_df['field_strength'].dropna()\nif len(fs) > 0:\n    fs_counts = fs.round(1).value_counts().sort_index()\n    axes[1].bar(fs_counts.index.astype(str), fs_counts.values, color='coral')\n    axes[1].set_xlabel(\"Field Strength (Tesla)\")\n    axes[1].set_ylabel(\"Count\")\n    axes[1].set_title(\"Magnetic Field Strength\", fontsize=13, fontweight='bold')\nelse:\n    axes[1].text(0.5, 0.5, \"No field strength data\", ha='center', va='center', transform=axes[1].transAxes)\n    axes[1].set_title(\"Magnetic Field Strength\", fontsize=13, fontweight='bold')\n\n# Bits stored\nbits_counts = dcm_df['bits_stored'].value_counts().sort_index()\naxes[2].bar(bits_counts.index.astype(str), bits_counts.values, color='mediumpurple')\naxes[2].set_xlabel(\"Bits Stored\")\naxes[2].set_ylabel(\"Count\")\naxes[2].set_title(\"Bits Per Pixel\", fontsize=13, fontweight='bold')\n\nplt.tight_layout()\nplt.show()\n\n# Laterality\nprint(f\"\\nLaterality distribution:\")\nprint(dcm_df['laterality'].value_counts().to_string())","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Pixel spacing vs slice thickness — important for resampling decisions\ndcm_df['plane'] = dcm_df['series_description'].apply(classify_plane)\n\nfig, ax = plt.subplots(figsize=(10, 7))\nfor plane, color in colors_plane.items():\n    mask = dcm_df['plane'] == plane\n    if mask.sum() > 0:\n        sub = dcm_df[mask].dropna(subset=['pixel_spacing_r', 'slice_thickness'])\n        ax.scatter(sub['pixel_spacing_r'], sub['slice_thickness'], \n                  alpha=0.4, s=20, c=color, label=f'{plane} ({mask.sum()})')\n\nax.set_xlabel(\"Pixel Spacing (mm)\", fontsize=12)\nax.set_ylabel(\"Slice Thickness (mm)\", fontsize=12)\nax.set_title(\"Pixel Spacing vs Slice Thickness by Anatomical Plane\", fontsize=14, fontweight='bold')\nax.legend(fontsize=10)\nax.set_xlim(0, None)\nax.set_ylim(0, None)\nplt.tight_layout()\nplt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n# 5. Image Visualization\n\nLet's look at actual knee MRI images across different planes, abnormalities, and windowing settings.","metadata":{}},{"cell_type":"code","source":"def load_dicom_image(dcm_path):\n    \"\"\"Load a DICOM file and return pixel array as float.\"\"\"\n    dcm = pydicom.dcmread(str(dcm_path))\n    img = dcm.pixel_array.astype(np.float32)\n    # Apply rescale if present\n    slope = getattr(dcm, 'RescaleSlope', 1)\n    intercept = getattr(dcm, 'RescaleIntercept', 0)\n    img = img * slope + intercept\n    return img, dcm\n\ndef apply_window(img, center, width):\n    \"\"\"Apply DICOM windowing.\"\"\"\n    lower = center - width / 2\n    upper = center + width / 2\n    img = np.clip(img, lower, upper)\n    img = (img - lower) / (upper - lower)\n    return img\n\ndef get_middle_slice(series_dir):\n    \"\"\"Get the middle slice from a series directory.\"\"\"\n    dcm_files = sorted(series_dir.glob(\"*.dcm\"))\n    if not dcm_files:\n        return None, None\n    mid_idx = len(dcm_files) // 2\n    return load_dicom_image(dcm_files[mid_idx])\n\nprint(\"Helper functions defined.\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Show sample images from different anatomical planes\n# Pick a study that has all 3 planes\nmerged = train_series_df.merge(train_df, on='StudyInstanceUID')\nmerged['plane'] = merged['Anatomical_Plane'].apply(classify_plane)\n\n# Find a study with all 3 main planes and some abnormalities\nstudy_plane_counts = merged.groupby('StudyInstanceUID')['plane'].apply(lambda x: set(x))\ngood_studies = study_plane_counts[study_plane_counts.apply(lambda x: {'Sagittal', 'Coronal', 'Axial'}.issubset(x))].index\n\n# Pick one with multiple abnormalities for visual interest\nabnormal_studies = labeled_df[labeled_df['n_abnormalities'] >= 3]\ncandidates = list(set(good_studies) & set(abnormal_studies['StudyInstanceUID']))\nif candidates:\n    demo_study_id = candidates[0]\nelse:\n    demo_study_id = good_studies[0] if len(good_studies) > 0 else train_studies[0].name\n\ndemo_study_dir = DATA_DIR / \"train_series\" / str(demo_study_id)\n\n# Get labels for this study\nstudy_labels = labeled_df[labeled_df['StudyInstanceUID'] == demo_study_id][target_cols].iloc[0]\npositive_labels = [col for col in target_cols if study_labels[col] == 1]\n\nprint(f\"Demo Study: {demo_study_id}\")\nprint(f\"Abnormalities: {positive_labels if positive_labels else 'Normal'}\")\nprint(f\"Number of series: {len(list(demo_study_dir.iterdir()))}\")\n\n# Show middle slice from each series in this study\nseries_dirs = sorted(demo_study_dir.iterdir())\nn_series = len(series_dirs)\nfig, axes = plt.subplots(1, min(n_series, 8), figsize=(4 * min(n_series, 8), 4))\nif n_series == 1:\n    axes = [axes]\n\nfor idx, series_dir in enumerate(series_dirs[:8]):\n    img, dcm = get_middle_slice(series_dir)\n    if img is None:\n        continue\n    \n    desc = getattr(dcm, 'SeriesDescription', 'Unknown')\n    plane = classify_plane(desc)\n    \n    # Normalize for display\n    p1, p99 = np.percentile(img, [1, 99])\n    img_display = np.clip((img - p1) / max(p99 - p1, 1), 0, 1)\n    \n    axes[idx].imshow(img_display, cmap='gray')\n    axes[idx].set_title(f\"{plane}\\n{desc[:30]}\", fontsize=9)\n    axes[idx].axis('off')\n\nfig.suptitle(f\"Study {demo_study_id} — All Series (middle slice)\\nLabels: {', '.join(positive_labels) if positive_labels else 'Normal'}\",\n             fontsize=12, fontweight='bold', y=1.05)\nplt.tight_layout()\nplt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Windowing Effects\nMRI windowing controls the brightness/contrast mapping. The same raw pixel data can reveal different structures depending on the window settings. This is critical for preprocessing — should we use DICOM window metadata or percentile-based normalization?","metadata":{}},{"cell_type":"code","source":"# Windowing comparison — same image with different windows\n# Pick a sagittal series for demonstration\nsag_series = merged[merged['plane'] == 'Sagittal'].iloc[0]\nsag_dir = DATA_DIR / \"train_series\" / str(sag_series['StudyInstanceUID'])\n\n# Find the sagittal series directory\nfor sd in sag_dir.iterdir():\n    if sd.is_dir():\n        dcm_files = sorted(sd.glob(\"*.dcm\"))\n        if dcm_files:\n            test_dcm = pydicom.dcmread(str(dcm_files[0]), stop_before_pixels=True)\n            desc = getattr(test_dcm, 'SeriesDescription', '')\n            if 'sag' in desc.lower():\n                img, dcm = load_dicom_image(dcm_files[len(dcm_files)//2])\n                break\n\nfig, axes = plt.subplots(1, 4, figsize=(20, 5))\n\n# Raw (no windowing)\naxes[0].imshow(img, cmap='gray')\naxes[0].set_title(f\"Raw\\nRange: [{img.min():.0f}, {img.max():.0f}]\", fontsize=10)\naxes[0].axis('off')\n\n# Percentile-based (1st-99th)\np1, p99 = np.percentile(img, [1, 99])\naxes[1].imshow(np.clip((img - p1) / max(p99 - p1, 1), 0, 1), cmap='gray')\naxes[1].set_title(f\"Percentile (1%-99%)\\n[{p1:.0f}, {p99:.0f}]\", fontsize=10)\naxes[1].axis('off')\n\n# DICOM window (if available)\nwc = dcm.WindowCenter if hasattr(dcm, 'WindowCenter') else img.mean()\nww = dcm.WindowWidth if hasattr(dcm, 'WindowWidth') else img.std() * 4\nif isinstance(wc, pydicom.multival.MultiValue): wc = float(wc[0])\nif isinstance(ww, pydicom.multival.MultiValue): ww = float(ww[0])\nwc, ww = float(wc), float(ww)\naxes[2].imshow(apply_window(img, wc, ww), cmap='gray')\naxes[2].set_title(f\"DICOM Window\\nC={wc:.0f}, W={ww:.0f}\", fontsize=10)\naxes[2].axis('off')\n\n# High contrast (narrow window)\naxes[3].imshow(apply_window(img, wc, ww * 0.5), cmap='gray')\naxes[3].set_title(f\"Narrow Window\\nC={wc:.0f}, W={ww*0.5:.0f}\", fontsize=10)\naxes[3].axis('off')\n\nfig.suptitle(\"Effect of Different Windowing on the Same MRI Slice\", fontsize=14, fontweight='bold', y=1.02)\nplt.tight_layout()\nplt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Normal vs. Abnormal — Can You Spot the Difference?\n\nLet's compare images from normal knees vs. knees with specific abnormalities. This helps build intuition for what the model needs to learn.","metadata":{}},{"cell_type":"code","source":"# Compare Normal vs Abnormal for key conditions\nconditions_to_show = ['ACL', 'Medial Meniscus', 'Effusion', 'Fracture']\n\nfig, axes = plt.subplots(len(conditions_to_show), 4, figsize=(16, 4 * len(conditions_to_show)))\n\nfor row_idx, condition in enumerate(conditions_to_show):\n    # Get 2 normal and 2 abnormal studies\n    normal_ids = labeled_df[labeled_df[condition] == 0]['StudyInstanceUID'].values\n    abnormal_ids = labeled_df[labeled_df[condition] == 1]['StudyInstanceUID'].values\n    \n    samples = []\n    for label, ids in [('Normal', normal_ids), ('Normal', normal_ids), \n                       (condition, abnormal_ids), (condition, abnormal_ids)]:\n        # Try to find a study with a sagittal series\n        found = False\n        for sid in np.random.RandomState(row_idx + len(samples)).permutation(ids)[:20]:\n            study_dir = DATA_DIR / \"train_series\" / str(sid)\n            if not study_dir.exists():\n                continue\n            for series_dir in study_dir.iterdir():\n                if not series_dir.is_dir():\n                    continue\n                dcm_files = sorted(series_dir.glob(\"*.dcm\"))\n                if not dcm_files:\n                    continue\n                test_dcm = pydicom.dcmread(str(dcm_files[0]), stop_before_pixels=True)\n                desc = getattr(test_dcm, 'SeriesDescription', '')\n                if 'sag' in desc.lower():\n                    img, _ = load_dicom_image(dcm_files[len(dcm_files)//2])\n                    samples.append((img, label, desc[:25]))\n                    found = True\n                    break\n            if found:\n                break\n        if not found:\n            samples.append((None, label, 'N/A'))\n    \n    for col_idx, (img, label, desc) in enumerate(samples):\n        ax = axes[row_idx, col_idx]\n        if img is not None:\n            p1, p99 = np.percentile(img, [1, 99])\n            img_display = np.clip((img - p1) / max(p99 - p1, 1), 0, 1)\n            ax.imshow(img_display, cmap='gray')\n        else:\n            ax.text(0.5, 0.5, 'No image', ha='center', va='center', transform=ax.transAxes)\n        \n        color = '#2ecc71' if label == 'Normal' else '#e74c3c'\n        ax.set_title(f\"{label}\", fontsize=10, color=color, fontweight='bold')\n        ax.axis('off')\n    \n    # Row label\n    axes[row_idx, 0].set_ylabel(condition, fontsize=12, fontweight='bold', rotation=0, \n                                 labelpad=60, va='center')\n\nfig.suptitle(\"Normal vs. Abnormal — Sagittal View Comparison\", fontsize=15, fontweight='bold', y=1.02)\nplt.tight_layout()\nplt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Slice-by-Slice Montage\nAn MRI series is a 3D volume. Let's visualize all slices in a single series to understand the spatial coverage.","metadata":{}},{"cell_type":"code","source":"# Montage — all slices from one series\n# Pick a series with a reasonable number of slices (20-40)\nmontage_study = demo_study_dir\nmontage_series = None\nfor sd in sorted(montage_study.iterdir()):\n    if sd.is_dir():\n        n_files = len(list(sd.glob(\"*.dcm\")))\n        if 15 <= n_files <= 40:\n            montage_series = sd\n            break\n\nif montage_series is None:\n    # Fallback: just use first series\n    montage_series = list(montage_study.iterdir())[0]\n\ndcm_files = sorted(montage_series.glob(\"*.dcm\"))\nn_slices = len(dcm_files)\nncols = 8\nnrows = (n_slices + ncols - 1) // ncols\n\n# Get series description\ntest_dcm = pydicom.dcmread(str(dcm_files[0]), stop_before_pixels=True)\ndesc = getattr(test_dcm, 'SeriesDescription', 'Unknown')\n\nfig, axes = plt.subplots(nrows, ncols, figsize=(ncols * 2.5, nrows * 2.5))\naxes = axes.ravel() if nrows > 1 else (axes if isinstance(axes, np.ndarray) else [axes])\n\nfor i in range(nrows * ncols):\n    if i < n_slices:\n        img, _ = load_dicom_image(dcm_files[i])\n        p1, p99 = np.percentile(img, [1, 99])\n        img_display = np.clip((img - p1) / max(p99 - p1, 1), 0, 1)\n        axes[i].imshow(img_display, cmap='gray')\n        axes[i].set_title(f\"Slice {i+1}\", fontsize=8)\n    axes[i].axis('off')\n\nfig.suptitle(f\"Complete Series Montage — {desc} ({n_slices} slices)\\nStudy: {demo_study_id}\",\n             fontsize=13, fontweight='bold', y=1.02)\nplt.tight_layout()\nplt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n# 6. Slice-Level Analysis\n\nHow many slices does each series have? This matters because models need to decide how many slices to sample from each series — too few and you miss pathology, too many and training slows down.","metadata":{}},{"cell_type":"code","source":"# Count slices per series across all training data\nslice_counts = []\nfor study_dir in train_studies:\n    if not study_dir.is_dir():\n        continue\n    study_id = study_dir.name\n    for series_dir in study_dir.iterdir():\n        if not series_dir.is_dir():\n            continue\n        n_dcm = len(list(series_dir.glob(\"*.dcm\")))\n        series_id = series_dir.name\n        \n        # Get plane from train_series_df\n        match = train_series_df[train_series_df['SeriesInstanceUID'] == series_id]\n        if len(match) > 0:\n            desc = match.iloc[0]['Anatomical_Plane']\n            plane = classify_plane(desc)\n        else:\n            plane = 'Unknown'\n        \n        slice_counts.append({\n            'StudyInstanceUID': study_id,\n            'SeriesInstanceUID': series_id,\n            'n_slices': n_dcm,\n            'plane': plane\n        })\n\nslices_df = pd.DataFrame(slice_counts)\nprint(f\"Total series analyzed: {len(slices_df):,}\")\nprint(f\"\\nSlice count statistics:\")\nprint(slices_df['n_slices'].describe().to_string())","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Slice count distributions by plane\nfig, axes = plt.subplots(1, 2, figsize=(16, 5))\n\n# Overall distribution\nax = axes[0]\nax.hist(slices_df['n_slices'], bins=60, color='steelblue', edgecolor='white', alpha=0.8)\nax.set_xlabel(\"Number of slices per series\", fontsize=12)\nax.set_ylabel(\"Count\", fontsize=12)\nax.set_title(\"Slices per Series — Overall Distribution\", fontsize=13, fontweight='bold')\nax.axvline(slices_df['n_slices'].median(), color='red', linestyle='--', \n           label=f\"Median: {slices_df['n_slices'].median():.0f}\")\nax.legend()\n\n# By plane (violin plot)\nax = axes[1]\nplanes_to_plot = ['Sagittal', 'Coronal', 'Axial']\nplot_data = slices_df[slices_df['plane'].isin(planes_to_plot)]\nif len(plot_data) > 0:\n    parts = ax.violinplot([plot_data[plot_data['plane'] == p]['n_slices'].values \n                          for p in planes_to_plot if (plot_data['plane'] == p).sum() > 0],\n                         showmeans=True, showmedians=True)\n    valid_planes = [p for p in planes_to_plot if (plot_data['plane'] == p).sum() > 0]\n    ax.set_xticks(range(1, len(valid_planes) + 1))\n    ax.set_xticklabels(valid_planes)\n    ax.set_ylabel(\"Number of slices\", fontsize=12)\n    ax.set_title(\"Slices per Series by Anatomical Plane\", fontsize=13, fontweight='bold')\n\nplt.tight_layout()\nplt.show()\n\n# Summary stats by plane\nprint(\"\\nSlices per series by plane:\")\nprint(slices_df.groupby('plane')['n_slices'].agg(['count', 'mean', 'median', 'min', 'max']).to_string())","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Pixel intensity distributions across planes\nfig, axes = plt.subplots(1, 3, figsize=(18, 5))\n\nrandom.seed(42)\nfor idx, (plane, ax) in enumerate(zip(['Sagittal', 'Coronal', 'Axial'], axes)):\n    plane_series = slices_df[slices_df['plane'] == plane]\n    if len(plane_series) == 0:\n        ax.text(0.5, 0.5, f'No {plane} series', ha='center', transform=ax.transAxes)\n        continue\n    \n    # Sample a few series and plot their intensity histograms\n    sample_ids = plane_series.sample(min(10, len(plane_series)), random_state=42)\n    for _, row in sample_ids.iterrows():\n        series_path = DATA_DIR / \"train_series\" / str(row['StudyInstanceUID']) / str(row['SeriesInstanceUID'])\n        dcm_files = sorted(series_path.glob(\"*.dcm\"))\n        if not dcm_files:\n            continue\n        img, _ = load_dicom_image(dcm_files[len(dcm_files)//2])\n        ax.hist(img.ravel(), bins=100, alpha=0.3, density=True)\n    \n    ax.set_xlabel(\"Pixel Intensity\")\n    ax.set_ylabel(\"Density\")\n    ax.set_title(f\"{plane} — Intensity Distributions\", fontsize=12, fontweight='bold')\n\nplt.suptitle(\"Pixel Intensity Distributions Across Planes (10 samples each)\",\n             fontsize=14, fontweight='bold', y=1.02)\nplt.tight_layout()\nplt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n# 7. Train vs Test Distribution\n\nA quick sanity check — do the train and test sets look similar in terms of series counts and metadata?","metadata":{}},{"cell_type":"code","source":"# Compare train vs test series metadata\ntest_series_df['plane'] = test_series_df['Anatomical_Plane'].apply(classify_plane)\n\nfig, axes = plt.subplots(1, 3, figsize=(18, 5))\n\n# Series per study\ntrain_sps = train_series_df.groupby('StudyInstanceUID').size()\ntest_sps = test_series_df.groupby('StudyInstanceUID').size()\n\nax = axes[0]\nax.hist(train_sps, bins=range(1, max(train_sps.max(), test_sps.max())+2), \n        alpha=0.6, color='steelblue', label=f'Train (n={len(train_sps)})', density=True, align='left')\nax.hist(test_sps, bins=range(1, max(train_sps.max(), test_sps.max())+2),\n        alpha=0.6, color='coral', label=f'Test (n={len(test_sps)})', density=True, align='left')\nax.set_xlabel(\"Series per study\")\nax.set_ylabel(\"Density\")\nax.set_title(\"Series per Study\", fontsize=13, fontweight='bold')\nax.legend()\n\n# Plane distribution comparison\nax = axes[1]\ntrain_planes = train_series_df['plane'].value_counts(normalize=True)\ntest_planes = test_series_df['plane'].value_counts(normalize=True)\nall_planes_union = sorted(set(train_planes.index) | set(test_planes.index))\nx = np.arange(len(all_planes_union))\nw = 0.35\nax.bar(x - w/2, [train_planes.get(p, 0) for p in all_planes_union], w, label='Train', color='steelblue')\nax.bar(x + w/2, [test_planes.get(p, 0) for p in all_planes_union], w, label='Test', color='coral')\nax.set_xticks(x)\nax.set_xticklabels(all_planes_union, rotation=30)\nax.set_ylabel(\"Fraction\")\nax.set_title(\"Plane Distribution\", fontsize=13, fontweight='bold')\nax.legend()\n\n# Series description overlap\nax = axes[2]\ntrain_descs = set(train_series_df['Anatomical_Plane'].dropna().unique())\ntest_descs = set(test_series_df['Anatomical_Plane'].dropna().unique())\noverlap = train_descs & test_descs\ntrain_only = train_descs - test_descs\ntest_only = test_descs - train_descs\n\nfrom matplotlib_venn import venn2\ntry:\n    venn2([train_descs, test_descs], set_labels=('Train', 'Test'), ax=ax)\n    ax.set_title(\"Series Description Overlap\", fontsize=13, fontweight='bold')\nexcept:\n    ax.bar(['Train only', 'Both', 'Test only'], [len(train_only), len(overlap), len(test_only)],\n           color=['steelblue', 'mediumpurple', 'coral'])\n    ax.set_title(\"Series Description Overlap\", fontsize=13, fontweight='bold')\n\nplt.tight_layout()\nplt.show()\n\nprint(f\"\\nTrain: {len(train_descs)} unique descriptions\")\nprint(f\"Test:  {len(test_descs)} unique descriptions\")\nprint(f\"Overlap: {len(overlap)}\")\nprint(f\"Test-only (unseen in train): {len(test_only)}\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n# 8. Full DICOM Header Reference\n\nFor reference, here's the complete DICOM header from one sample file. This is useful for identifying additional metadata fields you might want to leverage.","metadata":{}},{"cell_type":"code","source":"# Print full DICOM header from one sample file\nsample_study = train_studies[0]\nsample_series = list(sample_study.iterdir())[0]\nsample_dcm_file = list(sample_series.glob(\"*.dcm\"))[0]\n\ndcm = pydicom.dcmread(str(sample_dcm_file))\nprint(f\"DICOM Header — {sample_dcm_file.name}\")\nprint(\"=\" * 70)\nfor elem in dcm:\n    if elem.tag.group != 0x7FE0:  # Skip pixel data\n        print(f\"  {elem.tag} {elem.keyword:35s} {str(elem.value)[:80]}\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n# 9. Key Takeaways & Modeling Implications\n\n### Data Structure\n- Each **study** (patient) has multiple MRI **series** acquired in different anatomical planes\n- The 3 main planes are **Sagittal**, **Coronal**, and **Axial** — different abnormalities are best visible in different planes\n- Most studies have all 3 planes, but the exact number of series varies\n\n### Label Characteristics\n- **12 binary targets** — this is a **multi-label** classification problem\n- **Only 58 gold-labeled studies** out of 4,407 — heavily biased toward pathological cases (avg 4.14 abnormalities per knee)\n- **NLP label extraction** from radiology reports expands usable training data by ~75x\n- **Significant class imbalance** — rare conditions like Fracture and Contusion need special handling (weighted loss, oversampling)\n- **Strong co-occurrence patterns** — e.g., ACL tears often come with contusions and effusion. This suggests **multi-task learning** could help\n\n### Report Mining Insights\n- Reports are the primary label source — regex NLP achieves ~82% agreement with gold labels\n- Negation handling is critical (\"no ACL tear\" vs \"ACL tear\")\n- Some reports are multilingual — a dedicated multilingual pipeline could further improve coverage\n- NaN-filling strategy matters: treating unmentioned conditions as negative (0) is conservative but imperfect\n\n### Image Properties\n- **Variable image dimensions** — resizing/padding needed\n- **Variable pixel spacing** — consider resampling to uniform physical resolution\n- **Different scanners and field strengths** — domain adaptation/augmentation may help\n- **Percentile-based normalization** (1st-99th percentile) is generally safer than relying on DICOM window metadata\n\n### Modeling Strategy Hints\n1. **Multi-plane fusion**: The winning baseline (Pilkwang Kim, 524 votes) uses a \"slot\" system that maps each abnormality to its most relevant anatomical plane\n2. **Slice sampling**: Most series have 20-40 slices. The baseline samples 3 slices per plane (GROUP=3) and stacks as RGB\n3. **DINOv2 backbone**: Self-supervised ViT features transfer well to medical imaging\n4. **Anatomical priors**: ACL/Meniscus → Sagittal, OA → Coronal, Effusion → all planes\n5. **Label smoothing**: NLP-derived labels are noisy (~82% accuracy) — smooth toward 0.5 for NLP-only studies\n\n---\n**If this notebook was helpful, please give it an upvote!** It helps the community find useful resources.\n\nFeel free to fork this notebook and build on top of it. Happy kaggling!","metadata":{}},{"cell_type":"markdown","source":"---\n# 10. Baseline Model & Submission\n\nNow let's turn our EDA insights into a working model. Thanks to our **NLP label extraction** above, we're training on ~4,000+ studies instead of just 58.\n\n- **2.5D approach**: Stack 3 adjacent slices as RGB channels (the technique that wins RSNA competitions)\n- **ResNet18 backbone**: Small, fast, pretrained on ImageNet\n- **Multi-label BCE loss**: With class weights derived from expanded dataset prevalence\n- **Percentile normalization**: Based on our windowing analysis (safer than DICOM metadata)\n- **Sagittal-only**: We saw that sagittal is the most information-dense plane for most conditions\n\nThis baseline trains on Kaggle's free CPU and produces a submittable CSV.","metadata":{}},{"cell_type":"code","source":"import torch\nimport torch.nn as nn\nfrom torch.utils.data import Dataset, DataLoader\nfrom torchvision import models, transforms\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import roc_auc_score\nimport gc\n\n# Use CPU — expanded dataset is larger but still manageable on CPU\nDEVICE = torch.device('cpu')\nIMG_SIZE = 224\nNUM_SLICES = 3  # 2.5D: stack 3 adjacent slices as RGB\nBATCH_SIZE = 16\nEPOCHS = 8  # More epochs — we now have ~4,000+ studies instead of 58\nLR = 1e-4\n\nprint(f\"Device: {DEVICE}\")\nprint(f\"Image size: {IMG_SIZE}, Slices per input: {NUM_SLICES}\")\nprint(f\"Training with NLP-expanded labels: {len(labeled_df_expanded):,} studies (was {len(labeled_df)} gold-only)\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class KneeDataset(Dataset):\n    \"\"\"2.5D dataset: pick best sagittal series, sample 3 adjacent middle slices.\"\"\"\n    \n    def __init__(self, df, series_df, data_dir, target_cols, is_test=False):\n        self.df = df.reset_index(drop=True)\n        self.series_df = series_df\n        self.data_dir = Path(data_dir)\n        self.target_cols = target_cols\n        self.is_test = is_test\n        self.transform = transforms.Compose([\n            transforms.ToPILImage(),\n            transforms.Resize((IMG_SIZE, IMG_SIZE)),\n            transforms.ToTensor(),\n        ])\n        \n    def __len__(self):\n        return len(self.df)\n    \n    def _find_sagittal_series(self, study_id):\n        \"\"\"Find the best sagittal series for a study.\"\"\"\n        study_series = self.series_df[self.series_df['StudyInstanceUID'] == study_id]\n        sag = study_series[study_series['Anatomical_Plane'].str.lower().str.contains('sag', na=False)]\n        if len(sag) > 0:\n            return sag.iloc[0]['SeriesInstanceUID']\n        # Fallback: use first series\n        if len(study_series) > 0:\n            return study_series.iloc[0]['SeriesInstanceUID']\n        return None\n    \n    def _load_2_5d(self, study_id, series_id):\n        \"\"\"Load 3 adjacent middle slices as a 3-channel image.\"\"\"\n        folder = 'test_series' if self.is_test else 'train_series'\n        series_dir = self.data_dir / folder / str(study_id) / str(series_id)\n        \n        dcm_files = sorted(series_dir.glob(\"*.dcm\")) if series_dir.exists() else []\n        if not dcm_files:\n            return torch.zeros(3, IMG_SIZE, IMG_SIZE)\n        \n        mid = len(dcm_files) // 2\n        indices = [max(0, mid - 1), mid, min(len(dcm_files) - 1, mid + 1)]\n        \n        channels = []\n        for idx in indices:\n            dcm = pydicom.dcmread(str(dcm_files[idx]))\n            img = dcm.pixel_array.astype(np.float32)\n            slope = getattr(dcm, 'RescaleSlope', 1)\n            intercept = getattr(dcm, 'RescaleIntercept', 0)\n            img = img * slope + intercept\n            # Percentile normalization (from our EDA analysis)\n            p1, p99 = np.percentile(img, [1, 99])\n            img = np.clip((img - p1) / max(p99 - p1, 1), 0, 1)\n            img = (img * 255).astype(np.uint8)\n            channels.append(self.transform(img))\n        \n        return torch.cat(channels, dim=0)  # [3, H, W]\n    \n    def __getitem__(self, idx):\n        row = self.df.iloc[idx]\n        study_id = row['StudyInstanceUID']\n        \n        series_df = self.series_df[self.series_df['StudyInstanceUID'] == study_id]\n        series_id = self._find_sagittal_series(study_id)\n        \n        img = self._load_2_5d(study_id, series_id)\n        \n        if self.is_test:\n            return img, study_id\n        else:\n            labels = torch.tensor(row[self.target_cols].values.astype(np.float32))\n            return img, labels\n\nprint(\"KneeDataset defined.\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Model: ResNet18 with multi-label head\nclass KneeModel(nn.Module):\n    def __init__(self, num_classes=12):\n        super().__init__()\n        self.backbone = models.resnet18(weights='IMAGENET1K_V1')\n        in_features = self.backbone.fc.in_features\n        self.backbone.fc = nn.Sequential(\n            nn.Dropout(0.3),\n            nn.Linear(in_features, num_classes)\n        )\n    \n    def forward(self, x):\n        return self.backbone(x)\n\n# Compute class weights from EXPANDED dataset (not just 58 gold labels)\npos_counts = labeled_df_expanded[target_cols].sum()\nneg_counts = len(labeled_df_expanded) - pos_counts\npos_weight = (neg_counts / pos_counts.clip(lower=1)).values\npos_weight_tensor = torch.tensor(pos_weight, dtype=torch.float32).to(DEVICE)\n\nprint(f\"Model: ResNet18 → {len(target_cols)} outputs\")\nprint(f\"Class weights computed from {len(labeled_df_expanded):,} expanded studies\")\nprint(f\"\\nClass weights (neg/pos ratio):\")\nfor col, w in zip(target_cols, pos_weight):\n    print(f\"  {col:20s}: {w:.1f}\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Train/val split using EXPANDED labels (gold + NLP)\ntrain_data, val_data = train_test_split(labeled_df_expanded, test_size=0.10, random_state=42)\n\ntrain_dataset = KneeDataset(train_data, train_series_df, DATA_DIR, target_cols)\nval_dataset = KneeDataset(val_data, train_series_df, DATA_DIR, target_cols)\n\ntrain_loader = DataLoader(train_dataset, batch_size=BATCH_SIZE, shuffle=True, num_workers=0)\nval_loader = DataLoader(val_dataset, batch_size=BATCH_SIZE, shuffle=False, num_workers=0)\n\nmodel = KneeModel(num_classes=len(target_cols)).to(DEVICE)\ncriterion = nn.BCEWithLogitsLoss(pos_weight=pos_weight_tensor)\noptimizer = torch.optim.Adam(model.parameters(), lr=LR)\nscheduler = torch.optim.lr_scheduler.OneCycleLR(\n    optimizer, max_lr=LR, steps_per_epoch=len(train_loader), epochs=EPOCHS\n)\n\nprint(f\"Training on {len(train_data):,} studies, validating on {len(val_data):,} studies\")\nprint(f\"  (was {len(labeled_df)} total with gold-only — now {len(labeled_df_expanded):,} with NLP expansion)\")\nprint(f\"Train batches: {len(train_loader)}, Val batches: {len(val_loader)}\")\n\nimport sys\nbest_auc = 0\nfor epoch in range(EPOCHS):\n    # Train\n    model.train()\n    train_loss = 0\n    for imgs, labels in train_loader:\n        imgs, labels = imgs.to(DEVICE), labels.to(DEVICE)\n        optimizer.zero_grad()\n        logits = model(imgs)\n        loss = criterion(logits, labels)\n        loss.backward()\n        optimizer.step()\n        scheduler.step()\n        train_loss += loss.item()\n    \n    # Validate\n    model.eval()\n    val_preds, val_labels = [], []\n    val_loss = 0\n    with torch.no_grad():\n        for imgs, labels in val_loader:\n            imgs, labels = imgs.to(DEVICE), labels.to(DEVICE)\n            logits = model(imgs)\n            loss = criterion(logits, labels)\n            val_loss += loss.item()\n            val_preds.append(torch.sigmoid(logits).cpu().numpy())\n            val_labels.append(labels.cpu().numpy())\n    \n    val_preds = np.concatenate(val_preds)\n    val_labels = np.concatenate(val_labels)\n    \n    # Per-target AUC\n    aucs = []\n    for i, col in enumerate(target_cols):\n        if val_labels[:, i].sum() > 0 and val_labels[:, i].sum() < len(val_labels):\n            aucs.append(roc_auc_score(val_labels[:, i], val_preds[:, i]))\n    mean_auc = np.mean(aucs) if aucs else 0\n    \n    sys.stdout.flush()\n    print(f\"Epoch {epoch+1}/{EPOCHS} — \"\n          f\"Train Loss: {train_loss/len(train_loader):.4f}, \"\n          f\"Val Loss: {val_loss/len(val_loader):.4f}, \"\n          f\"Val AUC: {mean_auc:.4f}\")\n    \n    if mean_auc > best_auc:\n        best_auc = mean_auc\n        torch.save(model.state_dict(), '/kaggle/working/best_model.pth')\n\nprint(f\"\\nBest validation AUC: {best_auc:.4f}\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Generate submission\nmodel.load_state_dict(torch.load('/kaggle/working/best_model.pth', weights_only=True))\nmodel.eval()\n\ntest_dataset = KneeDataset(test_df, test_series_df, DATA_DIR, target_cols, is_test=True)\ntest_loader = DataLoader(test_dataset, batch_size=BATCH_SIZE, shuffle=False, num_workers=0)\n\nall_preds = {}\nwith torch.no_grad():\n    for imgs, study_ids in test_loader:\n        imgs = imgs.to(DEVICE)\n        logits = model(imgs)\n        probs = torch.sigmoid(logits).cpu().numpy()\n        for sid, pred in zip(study_ids, probs):\n            sid = sid if isinstance(sid, str) else str(sid.item())\n            all_preds[sid] = pred\n\n# Build submission in competition format\nsubmission = sample_sub.copy()\nfor i, col in enumerate(target_cols):\n    submission[col] = submission['StudyInstanceUID'].map(\n        lambda sid: all_preds.get(str(sid), np.full(len(target_cols), 0.5))[i]\n    )\n\nsubmission.to_csv('submission.csv', index=False)\nprint(f\"Submission saved: {submission.shape}\")\ndisplay(submission.head())\n\n# Cleanup\ndel model, train_dataset, val_dataset, test_dataset\ngc.collect()\npass  # no GPU cleanup needed","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n## What's Next? (Ideas to Improve This Baseline)\n\nThis baseline is deliberately minimal. Here's a ranked list of improvements by expected impact:\n\n1. **Use DINOv2 backbone** — the top-scoring approach (see Pilkwang Kim's 524-vote notebook)\n2. **Add coronal + axial planes** — our EDA showed different abnormalities are visible in different planes\n3. **Site-grouped cross-validation** — random folds leak scanner identity and inflate AUC by ~0.05\n4. **Improve NLP label quality** — multilingual patterns, sentence-level parsing, confidence scores\n5. **Label smoothing for NLP labels** — NLP labels are ~82% accurate; smooth toward 0.5 for noisy studies\n6. **More slices** — we only use 3 middle slices; try 5-7 or attention over all slices\n7. **Multi-plane attention head** — let the model learn which plane matters for each abnormality\n8. **Augmentation** — rigid transforms (rotation, translation, intensity scaling) — but NO horizontal flips (laterality matters!)\n\n---\n**If this notebook helped you, please upvote!** It helps the community find useful resources.\n\nFork this notebook to get started with your own experiments. Happy kaggling!","metadata":{}}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.10.0"}},"nbformat":4,"nbformat_minor":5}