{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.12.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceType":"competition","sourceId":24800,"datasetId":1042002,"databundleVersionId":1831594}],"dockerImageVersionId":31328,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# VinDr-CXR: Exploratory Data Analysis\nSystematic analysis of the VinDr-CXR chest X-ray dataset — class distribution, inter-rater agreement, soft labels, bounding box characteristics, and design implications for multi-label object detection with noisy annotations.","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport matplotlib.patches as mpatches\nimport seaborn as sns\nfrom pathlib import Path\nfrom collections import Counter\nfrom itertools import combinations\nimport warnings\n\nwarnings.filterwarnings('ignore')\nplt.rcParams.update({'font.size': 11, 'figure.dpi': 100})\nsns.set_style('whitegrid')\n\nINPUT_DIR = Path('/kaggle/input/competitions/vinbigdata-chest-xray-abnormalities-detection')\n\nCLASS_MAP = {\n    0: 'Aortic enlargement', 1: 'Atelectasis', 2: 'Calcification',\n    3: 'Cardiomegaly', 4: 'Consolidation', 5: 'ILD',\n    6: 'Infiltration', 7: 'Lung Opacity', 8: 'Nodule/Mass',\n    9: 'Other lesion', 10: 'Pleural effusion', 11: 'Pleural thickening',\n    12: 'Pneumothorax', 13: 'Pulmonary fibrosis', 14: 'No finding'\n}\n\nFOCAL_CLASSES = [1, 2, 5, 8, 9, 12]\nDIFFUSE_CLASSES = [0, 3, 4, 6, 7, 10, 11, 13]\n\nprint('Setup complete.')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:37:22.203078Z","iopub.execute_input":"2026-06-01T13:37:22.203426Z","iopub.status.idle":"2026-06-01T13:37:23.843101Z","shell.execute_reply.started":"2026-06-01T13:37:22.203393Z","shell.execute_reply":"2026-06-01T13:37:23.842214Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_df = pd.read_csv(INPUT_DIR / 'train.csv')\nprint(f'train.csv: {train_df.shape[0]:,} rows x {train_df.shape[1]} columns')\ntrain_df.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:37:23.844477Z","iopub.execute_input":"2026-06-01T13:37:23.844971Z","iopub.status.idle":"2026-06-01T13:37:24.020676Z","shell.execute_reply.started":"2026-06-01T13:37:23.844939Z","shell.execute_reply":"2026-06-01T13:37:24.019686Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 0.5 Data Preview: Normal vs. Abnormal X-rays\nUnderstanding the visual difference between a normal and abnormal chest X-ray is essential. A *Normal* image has no bounding boxes (class 14 = No finding), while an *Abnormal* image contains annotations from multiple radiologists — often disagreeing on both the presence and the class of findings. The overlay below shows each radiologist's bounding boxes in a different color, illustrating inter-rater uncertainty.","metadata":{}},{"cell_type":"code","source":"import pydicom\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\n\n\ndef read_xray(path, voi_lut=True, fix_monochrome=True):\n    dicom = pydicom.dcmread(path)\n    data = apply_voi_lut(dicom.pixel_array, dicom) if voi_lut else dicom.pixel_array\n    if fix_monochrome and dicom.PhotometricInterpretation == 'MONOCHROME1':\n        data = np.amax(data) - data\n    data = data - np.min(data)\n    data = data / np.max(data)\n    return (data * 255).astype(np.uint8)\n\n\n# Pick one Normal image (No finding) and one Abnormal image\nno_finding_imgs = train_df[train_df['class_id'] == 14]['image_id'].unique()\nfinding_imgs = train_df[train_df['class_id'] != 14]['image_id'].unique()\n\nnormal_id = no_finding_imgs[0]\nabnormal_id = finding_imgs[0]\n\nfig, axes = plt.subplots(1, 2, figsize=(14, 7))\n\n# Normal image\nimg_n = read_xray(str(INPUT_DIR / 'train' / f'{normal_id}.dicom'))\naxes[0].imshow(img_n, cmap='gray')\naxes[0].set_title('Normal (No Finding)', fontsize=13)\naxes[0].axis('off')\n\n# Abnormal image with bboxes from each radiologist\nimg_a = read_xray(str(INPUT_DIR / 'train' / f'{abnormal_id}.dicom'))\nabn_df = train_df[(train_df['image_id'] == abnormal_id) & (train_df['class_id'] != 14)]\nrad_ids = sorted(abn_df['rad_id'].unique())\nrad_colors = ['#e63946', '#2a9d8f', '#e9c46a', '#264653', '#f4a261', '#6a4c93']\n\naxes[1].imshow(img_a, cmap='gray')\nfor i, rid in enumerate(rad_ids[:6]):\n    rdf = abn_df[abn_df['rad_id'] == rid]\n    for _, row in rdf.iterrows():\n        x = row['x_min']\n        y = row['y_min']\n        w = row['x_max'] - row['x_min']\n        h = row['y_max'] - row['y_min']\n        rect = mpatches.Rectangle(\n            (x, y), w, h, linewidth=1.5,\n            edgecolor=rad_colors[i % len(rad_colors)],\n            facecolor='none'\n        )\n        axes[1].add_patch(rect)\n\nlegend_handles = [\n    mpatches.Patch(color=rad_colors[i % len(rad_colors)], label=f'Radiologist {rid}')\n    for i, rid in enumerate(rad_ids[:6])\n]\naxes[1].legend(handles=legend_handles, loc='best', fontsize=9)\naxes[1].set_title('Abnormal (multi-radiologist annotations)', fontsize=13)\naxes[1].axis('off')\n\nplt.tight_layout()\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:37:24.021904Z","iopub.execute_input":"2026-06-01T13:37:24.0223Z","iopub.status.idle":"2026-06-01T13:37:33.884671Z","shell.execute_reply.started":"2026-06-01T13:37:24.02227Z","shell.execute_reply":"2026-06-01T13:37:33.883559Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 1. Dataset Overview\nVerify structural integrity: column types, missing values, and basic sanity checks.","metadata":{}},{"cell_type":"code","source":"# 1. Dataset Overview\nprint('Shape:', train_df.shape)\nprint('\\nDtypes:')\nprint(train_df.dtypes)\nprint('\\nMissing values:')\nprint(train_df.isnull().sum())\nprint('\\nUnique images:', train_df['image_id'].nunique())\nprint('Unique radiologists:', train_df['rad_id'].nunique())\nprint('Unique classes:', train_df['class_id'].nunique())\n\n# Sanity: class_id range\nassert set(train_df['class_id'].unique()).issubset(set(range(15))), 'Unexpected class_id'\n# Sanity: No finding rows should have no bboxes\nnf = train_df[train_df['class_id'] == 14]\nprint(f'\\nNo-finding rows: {len(nf):,}')\nprint(f'  with non-null bbox: {nf[\"x_min\"].notna().sum()}')\n\n# Class name mapping\ntrain_df['class_name'] = train_df['class_id'].map(CLASS_MAP)\nprint('\\nValidation passed.')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:37:33.886859Z","iopub.execute_input":"2026-06-01T13:37:33.887195Z","iopub.status.idle":"2026-06-01T13:37:33.938184Z","shell.execute_reply.started":"2026-06-01T13:37:33.887165Z","shell.execute_reply":"2026-06-01T13:37:33.937181Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 2. Radiologist Analysis\nHow many radiologists annotated each image, and how are annotations distributed across them?","metadata":{}},{"cell_type":"code","source":"\n# 2. Radiologist Analysis\nrads = sorted(train_df['rad_id'].unique())\nR_POOL = len(rads)\nrads_per_img = train_df.groupby('image_id')['rad_id'].nunique()\n\nif rads_per_img.nunique() != 1:\n    print('WARNING: Radiologists per image is not constant. Using the mode for soft-label denominator.')\n\nR_PER_IMAGE = int(rads_per_img.mode().iloc[0])\nR = R_PER_IMAGE  # Backward-compatible alias used by later EDA cells.\n\nprint(f'Radiologist pool size: {R_POOL}')\nprint(f'Radiologists per image used for p_soft denominator (R): {R_PER_IMAGE}')\nprint(f'\\nRadiologists per image: mean={rads_per_img.mean():.2f}, '\n      f'min={rads_per_img.min()}, max={rads_per_img.max()}')\nprint(rads_per_img.value_counts().sort_index().to_frame('n_images'))\n\n# Workload distribution\nrad_workload = train_df.groupby('rad_id').size().sort_values(ascending=False)\nfig, ax = plt.subplots(figsize=(10, 4))\nrad_workload.plot(kind='bar', ax=ax, color='#457b9d')\nax.set_title('Annotations per Radiologist', fontsize=13)\nax.set_xlabel('Radiologist')\nax.set_ylabel('Number of annotations')\nplt.tight_layout()\nplt.show()\n\nprint(f'\\nWorkload imbalance ratio (max/min): {rad_workload.max() / rad_workload.min():.1f}x')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:37:33.939411Z","iopub.execute_input":"2026-06-01T13:37:33.939902Z","iopub.status.idle":"2026-06-01T13:37:34.313116Z","shell.execute_reply.started":"2026-06-01T13:37:33.939858Z","shell.execute_reply":"2026-06-01T13:37:34.312139Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 3. Class Distribution\nClass imbalance is a critical factor for loss function and sampling strategy design.","metadata":{}},{"cell_type":"code","source":"# 3. Class Distribution\nclass_counts = train_df.groupby(['class_id', 'class_name']).size().reset_index(name='count')\nclass_counts = class_counts.sort_values('count', ascending=False)\nn_images = train_df['image_id'].nunique()\nclass_counts['prevalence'] = class_counts['count'] / n_images\nprint(class_counts.to_string(index=False))\n\n# Plot\nfig, axes = plt.subplots(1, 2, figsize=(16, 5))\n\n# All classes\nsns.barplot(data=class_counts, y='class_name', x='count', ax=axes[0], palette='viridis')\naxes[0].set_title('Annotation Count per Class', fontsize=13)\naxes[0].set_xlabel('Count')\n\n# Prevalence (image-level)\nimg_class = train_df.drop_duplicates(subset=['image_id', 'class_id'])\nprev = img_class.groupby('class_name')['image_id'].count() / n_images\nprev = prev.sort_values(ascending=False)\nprev.plot(kind='barh', ax=axes[1], color='#e76f51')\naxes[1].set_title('Image-Level Prevalence', fontsize=13)\naxes[1].set_xlabel('Fraction of images')\n\nplt.tight_layout()\nplt.show()\n\nprint(f'\\nNo-finding prevalence: {prev.get(\"No finding\", 0):.2%}')\nprint(f'Most common abnormality: {prev.index[0] if prev.index[0] != \"No finding\" else prev.index[1]} ({prev.iloc[0] if prev.index[0] != \"No finding\" else prev.iloc[1]:.2%})')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:37:34.314255Z","iopub.execute_input":"2026-06-01T13:37:34.314714Z","iopub.status.idle":"2026-06-01T13:37:34.954282Z","shell.execute_reply.started":"2026-06-01T13:37:34.314683Z","shell.execute_reply":"2026-06-01T13:37:34.953265Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 4. Inter-Rater Agreement\nQuantify how often radiologists agree on the *presence* of any abnormality (Normal vs. Abnormal) and on specific class labels.","metadata":{}},{"cell_type":"code","source":"# 4a. Image-level consensus: Normal vs. Abnormal\nimg_abnormal = train_df.groupby('image_id').apply(\n    lambda g: (g['class_id'] != 14).any()\n).reset_index(name='any_abnormal')\n\n# For each image, count how many radiologists marked it as abnormal\nrad_abnormal = train_df.groupby(['image_id', 'rad_id']).apply(\n    lambda g: (g['class_id'] != 14).any()\n).reset_index(name='rad_sees_abnormal')\n\nconsensus = rad_abnormal.groupby('image_id')['rad_sees_abnormal'].agg(['sum', 'count'])\nconsensus.columns = ['n_abnormal_rads', 'n_total_rads']\nconsensus['n_normal_rads'] = consensus['n_total_rads'] - consensus['n_abnormal_rads']\n\n# Full consensus = all agree\nfull_consensus = (consensus['n_abnormal_rads'] == 0) | (consensus['n_abnormal_rads'] == consensus['n_total_rads'])\nprint('Image-level Normal/Abnormal consensus:')\nprint(f'  Full consensus: {full_consensus.sum()} / {len(consensus)} ({full_consensus.mean():.1%})')\nprint(f'  Disagreement: {(~full_consensus).sum()} images ({(~full_consensus).mean():.1%})')\n\n# Distribution\ndist = consensus['n_abnormal_rads'].value_counts().sort_index()\nprint(f'\\nDistribution of n_abnormal_rads per image:')\nprint(dist.to_string())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:37:34.955768Z","iopub.execute_input":"2026-06-01T13:37:34.956155Z","iopub.status.idle":"2026-06-01T13:37:42.945432Z","shell.execute_reply.started":"2026-06-01T13:37:34.956116Z","shell.execute_reply":"2026-06-01T13:37:42.944392Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 4b. Class-level agreement\nfindings_only = train_df[train_df['class_id'] != 14]\nimg_class_rad = findings_only.drop_duplicates(subset=['image_id', 'rad_id', 'class_id'])\n\n# For each (image, class): how many radiologists agree?\nclass_agree = img_class_rad.groupby(['image_id', 'class_id']).size().reset_index(name='n_rads_agree')\n\nagree_table = class_agree.groupby('class_id')['n_rads_agree'].agg(['mean', 'median', 'count'])\nagree_table.columns = ['mean_rads', 'median_rads', 'n_occurrences']\nagree_table['class_name'] = agree_table.index.map(CLASS_MAP)\nagree_table = agree_table[['class_name', 'n_occurrences', 'mean_rads', 'median_rads']].sort_values('mean_rads', ascending=False)\n\nprint('Class-level agreement (among radiologists who annotated the class):')\nprint(agree_table.to_string())\n\nprint(f'\\nClasses with mean agreement >= 2: {(agree_table[\"mean_rads\"] >= 2).sum()} / 14')\nprint(f'Classes with mean agreement < 1.5: {(agree_table[\"mean_rads\"] < 1.5).sum()} / 14')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:37:42.947024Z","iopub.execute_input":"2026-06-01T13:37:42.947501Z","iopub.status.idle":"2026-06-01T13:37:43.003167Z","shell.execute_reply.started":"2026-06-01T13:37:42.947457Z","shell.execute_reply":"2026-06-01T13:37:43.002128Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 5. Soft Label Analysis\nSoft labels capture annotation uncertainty: `p_soft = n_positive / R` represents the fraction of radiologists who labeled a class as present. This directly informs loss weighting and label smoothing strategies.","metadata":{}},{"cell_type":"code","source":"\n# 5a. Soft label calculation (vectorized)\nfindings_unique = train_df[train_df['class_id'] != 14].drop_duplicates(\n    subset=['image_id', 'rad_id', 'class_id']\n)\npositive_counts = findings_unique.groupby(['image_id', 'class_id']).size().reset_index(name='n_positive')\n\nall_img_ids = train_df['image_id'].unique()\nall_class_ids = list(range(14))\nfull_grid = pd.MultiIndex.from_product(\n    [all_img_ids, all_class_ids], names=['image_id', 'class_id']\n).to_frame(index=False)\n\nsoft_df = full_grid.merge(positive_counts, on=['image_id', 'class_id'], how='left')\nsoft_df['n_positive'] = soft_df['n_positive'].fillna(0).astype(int)\nsoft_df['p_soft'] = soft_df['n_positive'] / float(R_PER_IMAGE)\nsoft_df['class_name'] = soft_df['class_id'].map(CLASS_MAP)\n\nexpected_values = {i / R_PER_IMAGE for i in range(R_PER_IMAGE + 1)}\nobserved_values = set(np.round(soft_df['p_soft'].unique(), 10))\nif not observed_values.issubset({round(v, 10) for v in expected_values}):\n    raise RuntimeError(f'Unexpected p_soft values: {sorted(observed_values)}')\n\nprint(f'Soft label matrix: {soft_df.shape[0]:,} entries')\nprint(f'p_soft denominator R_PER_IMAGE: {R_PER_IMAGE}')\nprint(f'p_soft values: {sorted(soft_df[\"p_soft\"].unique())}')\nprint(f'p_soft > 0 entries: {(soft_df[\"p_soft\"] > 0).sum():,}')\nprint(f'p_soft == 1 entries: {(soft_df[\"p_soft\"] == 1.0).sum():,}')\nsoft_df.head(10)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:37:43.004828Z","iopub.execute_input":"2026-06-01T13:37:43.005252Z","iopub.status.idle":"2026-06-01T13:37:43.162034Z","shell.execute_reply.started":"2026-06-01T13:37:43.005207Z","shell.execute_reply":"2026-06-01T13:37:43.161128Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 5b. p_soft distribution per class (14 subplots)\nfig, axes = plt.subplots(2, 7, figsize=(18, 6))\naxes = axes.flatten()\n\nfor i in range(14):\n    ax = axes[i]\n    vals = soft_df[soft_df['class_id'] == i]['p_soft']\n    vals_nonzero = vals[vals > 0]\n    ax.hist(vals_nonzero, bins=np.arange(0.1, 1.1, 0.1), color='#457b9d', edgecolor='white')\n    ax.set_title(CLASS_MAP[i], fontsize=9)\n    ax.set_xlim(0, 1)\n    if i >= 7:\n        ax.set_xlabel('p_soft')\n    if i % 7 == 0:\n        ax.set_ylabel('Count')\n\nplt.suptitle('Soft Label Distribution per Class (p_soft > 0 only)', fontsize=13, y=1.02)\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:37:43.165424Z","iopub.execute_input":"2026-06-01T13:37:43.165763Z","iopub.status.idle":"2026-06-01T13:37:44.790852Z","shell.execute_reply.started":"2026-06-01T13:37:43.165733Z","shell.execute_reply":"2026-06-01T13:37:44.789676Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 5c. Summary table with disagreement rates\nsoft_summary = []\nfor cid in range(14):\n    vals = soft_df[soft_df['class_id'] == cid]['p_soft']\n    nonzero = vals[vals > 0]\n    soft_summary.append({\n        'class_name': CLASS_MAP[cid],\n        'n_images_positive': len(nonzero),\n        'p_soft_mean': nonzero.mean() if len(nonzero) > 0 else 0,\n        'p_soft_median': nonzero.median() if len(nonzero) > 0 else 0,\n        'full_agree_pct': (nonzero == 1.0).sum() / max(len(nonzero), 1),\n        'disagree_pct': ((nonzero > 0) & (nonzero < 1.0)).sum() / max(len(nonzero), 1)\n    })\n\nsoft_summary_df = pd.DataFrame(soft_summary).sort_values('disagree_pct', ascending=False)\nprint(soft_summary_df.to_string(index=False))\n\nhigh_disagree = soft_summary_df[soft_summary_df['disagree_pct'] > 0.5]\nprint(f'\\nClasses with >50% disagreement: {len(high_disagree)} / 14')\nprint('These classes benefit most from soft labels / label smoothing.')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:37:44.791943Z","iopub.execute_input":"2026-06-01T13:37:44.792224Z","iopub.status.idle":"2026-06-01T13:37:44.851584Z","shell.execute_reply.started":"2026-06-01T13:37:44.792197Z","shell.execute_reply":"2026-06-01T13:37:44.850434Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 6. Resolution & DICOM Metadata\nImage resolution affects detection anchor design and resizing strategy. DICOM metadata (Window Center/Width, RescaleSlope, PhotometricInterpretation) informs preprocessing decisions.","metadata":{}},{"cell_type":"code","source":"\n# 6a. DICOM metadata extraction\nimport os\nimport pydicom\n\ntrain_dir = INPUT_DIR / 'train'\nall_metadata_ids = train_df['image_id'].unique()\nmetadata_sample_size = int(os.environ.get('EDA_DICOM_METADATA_SAMPLE_SIZE', '0') or '0')\nif metadata_sample_size > 0:\n    sample_ids = all_metadata_ids[:metadata_sample_size]\n    print(f'Extracting DICOM metadata sample: {metadata_sample_size:,}/{len(all_metadata_ids):,} images')\nelse:\n    sample_ids = all_metadata_ids\n    print(f'Extracting DICOM metadata for all {len(sample_ids):,} images')\n\nmeta_rows = []\nmetadata_errors = []\nfor img_id in sample_ids:\n    dcm_path = train_dir / f'{img_id}.dicom'\n    if not os.path.exists(dcm_path):\n        metadata_errors.append({'image_id': img_id, 'error': 'missing_dicom'})\n        continue\n    try:\n        ds = pydicom.dcmread(str(dcm_path), stop_before_pixels=True)\n        meta_rows.append({\n            'image_id': img_id,\n            'Rows': int(ds.Rows),\n            'Columns': int(ds.Columns),\n            'BitsAllocated': int(ds.BitsAllocated) if hasattr(ds, 'BitsAllocated') else None,\n            'BitsStored': int(ds.BitsStored) if hasattr(ds, 'BitsStored') else None,\n            'PixelSpacing': float(ds.PixelSpacing[0]) if hasattr(ds, 'PixelSpacing') else None,\n            'RescaleSlope': float(ds.RescaleSlope) if hasattr(ds, 'RescaleSlope') else None,\n            'RescaleIntercept': float(ds.RescaleIntercept) if hasattr(ds, 'RescaleIntercept') else None,\n            'PhotometricInterpretation': str(ds.PhotometricInterpretation) if hasattr(ds, 'PhotometricInterpretation') else None,\n            'WindowCenter': str(ds.WindowCenter) if hasattr(ds, 'WindowCenter') else None,\n            'WindowWidth': str(ds.WindowWidth) if hasattr(ds, 'WindowWidth') else None,\n        })\n    except Exception as exc:\n        metadata_errors.append({'image_id': img_id, 'error': type(exc).__name__})\n\nmeta_df = pd.DataFrame(meta_rows)\nmetadata_error_df = pd.DataFrame(metadata_errors)\nprint(f'Extracted metadata for {len(meta_df):,}/{len(sample_ids):,} requested images')\nif len(metadata_error_df):\n    print('Metadata errors:')\n    print(metadata_error_df['error'].value_counts().to_string())\nmeta_df.head()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:37:44.852928Z","iopub.execute_input":"2026-06-01T13:37:44.853512Z","iopub.status.idle":"2026-06-01T13:41:02.285344Z","shell.execute_reply.started":"2026-06-01T13:37:44.853468Z","shell.execute_reply":"2026-06-01T13:41:02.284119Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 6b. Resolution statistics\nprint('Resolution statistics:')\nprint(meta_df[['Rows', 'Columns']].describe())\n\nfig, axes = plt.subplots(1, 2, figsize=(12, 4))\nmeta_df['Rows'].hist(ax=axes[0], bins=30, color='#457b9d')\naxes[0].set_title('Height Distribution', fontsize=13)\naxes[0].set_xlabel('Pixels')\n\nmeta_df['Columns'].hist(ax=axes[1], bins=30, color='#e76f51')\naxes[1].set_title('Width Distribution', fontsize=13)\naxes[1].set_xlabel('Pixels')\n\nplt.tight_layout()\nplt.show()\n\nprint(f'\\nMedian resolution: {meta_df[\"Rows\"].median():.0f} x {meta_df[\"Columns\"].median():.0f}')\nprint(f'Common aspect ratio: {(meta_df[\"Columns\"] / meta_df[\"Rows\"]).median():.2f}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:41:02.28661Z","iopub.execute_input":"2026-06-01T13:41:02.286904Z","iopub.status.idle":"2026-06-01T13:41:02.681439Z","shell.execute_reply.started":"2026-06-01T13:41:02.286875Z","shell.execute_reply":"2026-06-01T13:41:02.680491Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 6c. Preprocessing metadata\nprint('PhotometricInterpretation:')\nprint(meta_df['PhotometricInterpretation'].value_counts())\n\nprint(f'\\nRescaleSlope range: {meta_df[\"RescaleSlope\"].min()} - {meta_df[\"RescaleSlope\"].max()}')\nprint(f'PixelSpacing range: {meta_df[\"PixelSpacing\"].min():.2f} - {meta_df[\"PixelSpacing\"].max():.2f} mm')\n\n# Check for WindowCenter / WindowWidth (VOI LUT)\nwc_vals, ww_vals = [], []\nfor img_id in sample_ids[:200]:\n    dcm_path = train_dir / f'{img_id}.dicom'\n    if not os.path.exists(dcm_path):\n        continue\n    ds = pydicom.dcmread(str(dcm_path), stop_before_pixels=True)\n    if hasattr(ds, 'WindowCenter'):\n        wc = ds.WindowCenter\n        wc_vals.append(float(wc) if not isinstance(wc, pydicom.valuerep.DSfloat) else float(wc))\n    if hasattr(ds, 'WindowWidth'):\n        ww = ds.WindowWidth\n        ww_vals.append(float(ww) if not isinstance(ww, pydicom.valuerep.DSfloat) else float(ww))\n\nif wc_vals:\n    print(f'\\nWindowCenter: mean={np.mean(wc_vals):.0f}, range=[{np.min(wc_vals):.0f}, {np.max(wc_vals):.0f}]')\nif ww_vals:\n    print(f'WindowWidth: mean={np.mean(ww_vals):.0f}, range=[{np.min(ww_vals):.0f}, {np.max(ww_vals):.0f}]')\n\nprint('\\nRecommendation: Use apply_voi_lut=True and fix_monochrome=True in read_xray().')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:41:02.682707Z","iopub.execute_input":"2026-06-01T13:41:02.683094Z","iopub.status.idle":"2026-06-01T13:41:03.186682Z","shell.execute_reply.started":"2026-06-01T13:41:02.683052Z","shell.execute_reply":"2026-06-01T13:41:03.185718Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 7. Bounding Box Analysis\nBounding box size distributions inform anchor design and feature map stride choices. Small lesions require finer feature maps; diffuse findings may need larger receptive fields.","metadata":{}},{"cell_type":"code","source":"\n# 7a. Box statistics per class and area analysis\nbbox_df = train_df[train_df['class_id'] != 14].dropna(subset=['x_min', 'y_min', 'x_max', 'y_max']).copy()\nbbox_df['width'] = bbox_df['x_max'] - bbox_df['x_min']\nbbox_df['height'] = bbox_df['y_max'] - bbox_df['y_min']\nbbox_df['area'] = bbox_df['width'] * bbox_df['height']\n\nif 'meta_df' in globals() and len(meta_df) > 0 and {'Rows', 'Columns'}.issubset(meta_df.columns):\n    bbox_df = bbox_df.merge(meta_df[['image_id', 'Rows', 'Columns']], on='image_id', how='left')\n    bbox_df['area_ratio'] = bbox_df['area'] / (bbox_df['Rows'] * bbox_df['Columns'])\n    bbox_df['width_at_1024'] = bbox_df['width'] * (1024.0 / bbox_df['Columns'])\n    bbox_df['height_at_1024'] = bbox_df['height'] * (1024.0 / bbox_df['Rows'])\nelse:\n    bbox_df['area_ratio'] = np.nan\n    bbox_df['width_at_1024'] = np.nan\n    bbox_df['height_at_1024'] = np.nan\n    print('WARNING: DICOM resolution metadata is unavailable. Normalized box-size metrics will be NaN.')\n\nbbox_stats = []\nfor cid in range(14):\n    cls_data = bbox_df[bbox_df['class_id'] == cid]\n    if len(cls_data) == 0:\n        continue\n    median_min_side_cells_s16 = min(\n        cls_data['width_at_1024'].median(skipna=True),\n        cls_data['height_at_1024'].median(skipna=True),\n    ) / 16.0\n    bbox_stats.append({\n        'class_id': cid,\n        'class_name': CLASS_MAP[cid],\n        'n_boxes': int(len(cls_data)),\n        'n_images': int(cls_data['image_id'].nunique()),\n        'area_median': float(cls_data['area'].median()),\n        'area_ratio_median': float(cls_data['area_ratio'].median(skipna=True)) if cls_data['area_ratio'].notna().any() else np.nan,\n        'area_ratio_p10': float(cls_data['area_ratio'].quantile(0.10)) if cls_data['area_ratio'].notna().any() else np.nan,\n        'width_median': float(cls_data['width'].median()),\n        'height_median': float(cls_data['height'].median()),\n        'width_at_1024_median': float(cls_data['width_at_1024'].median(skipna=True)) if cls_data['width_at_1024'].notna().any() else np.nan,\n        'height_at_1024_median': float(cls_data['height_at_1024'].median(skipna=True)) if cls_data['height_at_1024'].notna().any() else np.nan,\n        'median_min_side_cells_s16': float(median_min_side_cells_s16) if pd.notna(median_min_side_cells_s16) else np.nan,\n        'group': 'focal' if cid in FOCAL_CLASSES else 'diffuse',\n    })\n\nbbox_stats_df = pd.DataFrame(bbox_stats).sort_values('area_ratio_median', ascending=False, na_position='last')\nprint(bbox_stats_df.to_string(index=False))\n\nfig, ax = plt.subplots(figsize=(14, 5))\nplot_df = bbox_stats_df.sort_values('area_ratio_median', ascending=True, na_position='last')\ncolors = ['#e76f51' if cid in FOCAL_CLASSES else '#457b9d' for cid in plot_df['class_id']]\nax.barh(plot_df['class_name'], plot_df['area_ratio_median'] * 100.0, color=colors, alpha=0.75)\nax.set_title('Median Box Area Ratio — Focal=orange, Diffuse=blue', fontsize=13)\nax.set_xlabel('Median bbox area / image area (%)')\nplt.tight_layout()\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:41:03.187938Z","iopub.execute_input":"2026-06-01T13:41:03.188253Z","iopub.status.idle":"2026-06-01T13:41:03.548506Z","shell.execute_reply.started":"2026-06-01T13:41:03.188216Z","shell.execute_reply":"2026-06-01T13:41:03.547606Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# 7b. Small lesion analysis + feature map cell calculation\nINPUT_SIZE = 1024\nSTRIDES = [8, 16, 32, 64]\n\nsmall_stats = []\nfor cid in range(14):\n    class_boxes = bbox_df[bbox_df['class_id'] == cid]\n    if len(class_boxes) == 0:\n        continue\n    if class_boxes['width_at_1024'].notna().any() and class_boxes['height_at_1024'].notna().any():\n        scaled_width = class_boxes['width_at_1024']\n        scaled_height = class_boxes['height_at_1024']\n    else:\n        # Fallback only for environments where DICOM metadata was not available.\n        scaled_width = class_boxes['width'] * (INPUT_SIZE / 3000.0)\n        scaled_height = class_boxes['height'] * (INPUT_SIZE / 3000.0)\n    for stride in STRIDES:\n        small_stats.append({\n            'class_id': cid,\n            'class': CLASS_MAP[cid],\n            'stride': stride,\n            'pct_width_lt_cell': float((scaled_width < stride).mean()),\n            'pct_height_lt_cell': float((scaled_height < stride).mean()),\n            'p10_min_side_cells': float(min(scaled_width.quantile(0.10), scaled_height.quantile(0.10)) / stride),\n            'median_min_side_cells': float(min(scaled_width.median(), scaled_height.median()) / stride),\n        })\n\nsmall_df = pd.DataFrame(small_stats)\nprint('Fraction of boxes smaller than feature map cell at each stride:')\npivot = small_df.pivot(index='class', columns='stride', values='pct_width_lt_cell')\npivot = pivot.sort_values(pivot.columns[0], ascending=False)\nprint(pivot.to_string(float_format='%.2f'))\n\nstride16 = small_df[small_df['stride'] == 16].sort_values('median_min_side_cells')\nprint('\\nMedian min-side cells at stride 16 after resizing to 1024:')\nprint(stride16[['class', 'p10_min_side_cells', 'median_min_side_cells']].to_string(index=False, float_format='%.2f'))\n\nprint('\\nInterpretation: detector routing should penalize classes whose boxes remain below a few feature-map cells after resize.')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:41:03.550054Z","iopub.execute_input":"2026-06-01T13:41:03.550509Z","iopub.status.idle":"2026-06-01T13:41:03.688055Z","shell.execute_reply.started":"2026-06-01T13:41:03.55047Z","shell.execute_reply":"2026-06-01T13:41:03.687207Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 8. IoU Analysis\nInter-rater IoU quantifies spatial agreement. Low IoU classes may need relaxed matching thresholds (`theta_IoU`) during training to avoid penalizing reasonable predictions.","metadata":{}},{"cell_type":"code","source":"\n# 8a. Per-class IoU calculation\n\ndef compute_iou(box1, box2):\n    x1 = max(box1[0], box2[0])\n    y1 = max(box1[1], box2[1])\n    x2 = min(box1[2], box2[2])\n    y2 = min(box1[3], box2[3])\n    inter = max(0, x2 - x1) * max(0, y2 - y1)\n    a1 = max(0, box1[2] - box1[0]) * max(0, box1[3] - box1[1])\n    a2 = max(0, box2[2] - box2[0]) * max(0, box2[3] - box2[1])\n    union = a1 + a2 - inter\n    return inter / union if union > 0 else 0.0\n\n# Symmetric best-match IoU per image/class/radiologist-pair.\n# This avoids the older all-pairs product that over-penalized cases with multiple boxes.\niou_records = []\nfor cid in range(14):\n    class_bboxes = bbox_df[bbox_df['class_id'] == cid]\n    for img_id, grp in class_bboxes.groupby('image_id'):\n        rads_in_img = sorted(grp['rad_id'].unique())\n        if len(rads_in_img) < 2:\n            continue\n        rad_boxes = {\n            rid: list(zip(\n                grp[grp['rad_id'] == rid]['x_min'],\n                grp[grp['rad_id'] == rid]['y_min'],\n                grp[grp['rad_id'] == rid]['x_max'],\n                grp[grp['rad_id'] == rid]['y_max'],\n            ))\n            for rid in rads_in_img\n        }\n        for r1, r2 in combinations(rads_in_img, 2):\n            boxes1 = rad_boxes[r1]\n            boxes2 = rad_boxes[r2]\n            if not boxes1 or not boxes2:\n                continue\n            for b1 in boxes1:\n                iou_records.append({'class_id': cid, 'class_name': CLASS_MAP[cid], 'iou': max(compute_iou(b1, b2) for b2 in boxes2)})\n            for b2 in boxes2:\n                iou_records.append({'class_id': cid, 'class_name': CLASS_MAP[cid], 'iou': max(compute_iou(b2, b1) for b1 in boxes1)})\n\niou_df = pd.DataFrame(iou_records)\nprint(f'Total symmetric best-match IoU computations: {len(iou_df):,}')\niou_df.head()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:41:03.689285Z","iopub.execute_input":"2026-06-01T13:41:03.689555Z","iopub.status.idle":"2026-06-01T13:41:36.337651Z","shell.execute_reply.started":"2026-06-01T13:41:03.689527Z","shell.execute_reply":"2026-06-01T13:41:36.336648Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# 8b. IoU statistics table + theta_IoU recommendations\nif iou_df.empty:\n    raise RuntimeError('No pairwise IoU records were produced. Check bbox_df and radiologist annotations.')\n\niou_stats = iou_df.groupby('class_id')['iou'].agg(\n    n_comparisons='size',\n    mean='mean',\n    median='median',\n    std='std',\n    p25=lambda x: x.quantile(0.25),\n    p40=lambda x: x.quantile(0.40),\n    p50=lambda x: x.quantile(0.50),\n    p75=lambda x: x.quantile(0.75),\n)\niou_stats['class_name'] = iou_stats.index.map(CLASS_MAP)\niou_stats['pct_iou_0'] = iou_df.groupby('class_id')['iou'].apply(lambda x: float((x == 0).mean() * 100.0))\niou_stats['pct_iou_lt_03'] = iou_df.groupby('class_id')['iou'].apply(lambda x: float((x < 0.3).mean() * 100.0))\niou_stats['theta_iou_rec'] = iou_stats['p40'].clip(lower=0.15, upper=0.50).round(2)\niou_stats = iou_stats[['class_name', 'n_comparisons', 'mean', 'median', 'p25', 'p40', 'p75', 'pct_iou_0', 'pct_iou_lt_03', 'theta_iou_rec']]\n\ndisplay_iou = iou_stats.sort_values('median', ascending=False)\nprint('Per-class symmetric best-match IoU statistics:')\nprint(display_iou.to_string())\n\nprint(f'\\nOverall mean IoU: {iou_df[\"iou\"].mean():.3f}')\nprint(f'Overall median IoU: {iou_df[\"iou\"].median():.3f}')\nlow_iou = iou_stats[iou_stats['median'] < 0.30]\nprint(f'Classes with median IoU < 0.30: {len(low_iou)} (weak/conditional detector candidates)')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:41:36.338874Z","iopub.execute_input":"2026-06-01T13:41:36.339225Z","iopub.status.idle":"2026-06-01T13:41:36.412371Z","shell.execute_reply.started":"2026-06-01T13:41:36.339176Z","shell.execute_reply":"2026-06-01T13:41:36.411263Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 8c. IoU distribution boxplot\nfig, ax = plt.subplots(figsize=(14, 5))\norder = iou_stats.sort_values('mean').index.tolist()\nplot_data = iou_df[iou_df['class_id'].isin(order)]\nplot_data['class_name'] = plot_data['class_id'].map(CLASS_MAP)\n\nsns.boxplot(data=plot_data, x='class_name', y='iou', order=[CLASS_MAP[c] for c in order],\n            ax=ax, palette='Set2', fliersize=1)\nax.set_title('Inter-Rater IoU Distribution per Class', fontsize=13)\nax.set_xlabel('')\nax.set_ylabel('IoU')\nax.axhline(y=0.5, color='red', linestyle='--', alpha=0.5, label='IoU=0.5')\nax.axhline(y=0.25, color='orange', linestyle='--', alpha=0.5, label='IoU=0.25')\nax.legend(loc='best')\nplt.xticks(rotation=45, ha='right')\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:41:36.413582Z","iopub.execute_input":"2026-06-01T13:41:36.413937Z","iopub.status.idle":"2026-06-01T13:41:37.116744Z","shell.execute_reply.started":"2026-06-01T13:41:36.413897Z","shell.execute_reply":"2026-06-01T13:41:37.115718Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 9. WBF Analysis\nWeighted Boxes Fusion (WBF) aggregates multi-radiologist annotations into consensus boxes. The IoU threshold parameter controls the trade-off between retention and coverage.","metadata":{}},{"cell_type":"code","source":"!pip install ensemble-boxes","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:41:37.117891Z","iopub.execute_input":"2026-06-01T13:41:37.118254Z","iopub.status.idle":"2026-06-01T13:41:43.08345Z","shell.execute_reply.started":"2026-06-01T13:41:37.118215Z","shell.execute_reply":"2026-06-01T13:41:43.082276Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# 9a. WBF fusion\ntry:\n    from ensemble_boxes import weighted_boxes_fusion\nexcept ImportError as exc:\n    raise ImportError('Install ensemble-boxes before running WBF analysis: pip install ensemble-boxes') from exc\n\nwbf_source_df = bbox_df.copy()\nWBF_IOU_THR = 0.5\nWBF_SKIP_BOX_THR = 0.0\n\ndef _image_shape_for_wbf(group):\n    if {'Rows', 'Columns'}.issubset(group.columns) and group['Rows'].notna().any() and group['Columns'].notna().any():\n        return float(group['Rows'].dropna().iloc[0]), float(group['Columns'].dropna().iloc[0])\n    return 3000.0, 3000.0\n\nwbf_records = []\nprint('Running WBF by (image_id, class_id)...')\nfor (img_id, cid), group in wbf_source_df.groupby(['image_id', 'class_id']):\n    img_h, img_w = _image_shape_for_wbf(group)\n    boxes_list, scores_list, labels_list = [], [], []\n    for rad_id in sorted(group['rad_id'].unique()):\n        rad_df = group[group['rad_id'] == rad_id]\n        boxes = rad_df[['x_min', 'y_min', 'x_max', 'y_max']].astype(float).copy()\n        boxes['x_min'] = (boxes['x_min'] / img_w).clip(0.0, 1.0)\n        boxes['x_max'] = (boxes['x_max'] / img_w).clip(0.0, 1.0)\n        boxes['y_min'] = (boxes['y_min'] / img_h).clip(0.0, 1.0)\n        boxes['y_max'] = (boxes['y_max'] / img_h).clip(0.0, 1.0)\n        valid = boxes[(boxes['x_max'] > boxes['x_min']) & (boxes['y_max'] > boxes['y_min'])]\n        boxes_list.append(valid[['x_min', 'y_min', 'x_max', 'y_max']].values.tolist())\n        scores_list.append([1.0] * len(valid))\n        labels_list.append([int(cid)] * len(valid))\n    if not any(len(items) for items in boxes_list):\n        continue\n    fused_boxes, fused_scores, fused_labels = weighted_boxes_fusion(\n        boxes_list,\n        scores_list,\n        labels_list,\n        iou_thr=WBF_IOU_THR,\n        skip_box_thr=WBF_SKIP_BOX_THR,\n    )\n    for box, score, label in zip(fused_boxes, fused_scores, fused_labels):\n        wbf_records.append({\n            'image_id': img_id,\n            'class_id': int(label),\n            'class_name': CLASS_MAP[int(label)],\n            'wbf_score': float(score),\n            'x1_norm': float(box[0]),\n            'y1_norm': float(box[1]),\n            'x2_norm': float(box[2]),\n            'y2_norm': float(box[3]),\n            'source_box_count': int(len(group)),\n        })\n\nwbf_df_results = pd.DataFrame(wbf_records)\nprint(f'WBF complete: {len(wbf_df_results):,} fused boxes from {len(wbf_source_df):,} source boxes')\nwbf_df_results.head()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:41:43.085358Z","iopub.execute_input":"2026-06-01T13:41:43.085769Z","iopub.status.idle":"2026-06-01T13:44:53.776746Z","shell.execute_reply.started":"2026-06-01T13:41:43.085731Z","shell.execute_reply":"2026-06-01T13:44:53.774404Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# 9b. Retention / Coverage table\nraw_stats = bbox_df.groupby('class_id').agg(\n    raw_boxes=('class_id', 'size'),\n    raw_images=('image_id', 'nunique'),\n).reset_index()\nraw_stats['class_name'] = raw_stats['class_id'].map(CLASS_MAP)\n\nthresholds = {\n    f'>=1/{R_PER_IMAGE}': (1.0 / R_PER_IMAGE) * 0.95,\n    f'>=2/{R_PER_IMAGE}': (2.0 / R_PER_IMAGE) * 0.95,\n    f'>=3/{R_PER_IMAGE}': (3.0 / R_PER_IMAGE) * 0.95,\n}\n\nwbf_analysis = []\nfor _, row in raw_stats.iterrows():\n    cid = int(row['class_id'])\n    cls_wbf = wbf_df_results[wbf_df_results['class_id'] == cid] if len(wbf_df_results) else pd.DataFrame()\n    result = {\n        'class_id': cid,\n        'class_name': row['class_name'],\n        'raw_boxes': int(row['raw_boxes']),\n        'raw_images': int(row['raw_images']),\n        'group': 'focal' if cid in FOCAL_CLASSES else 'diffuse',\n    }\n    for name, threshold in thresholds.items():\n        valid = cls_wbf[cls_wbf['wbf_score'] >= threshold] if len(cls_wbf) else cls_wbf\n        n_boxes = int(len(valid))\n        n_images = int(valid['image_id'].nunique()) if len(valid) else 0\n        result[f'wbf_boxes_{name}'] = n_boxes\n        result[f'wbf_retention_{name}'] = n_boxes / max(int(row['raw_boxes']), 1)\n        result[f'wbf_cover_{name}'] = n_images / max(int(row['raw_images']), 1)\n    wbf_analysis.append(result)\n\nwbf_di_df = pd.DataFrame(wbf_analysis)\nwbf_analysis_df = wbf_di_df.sort_values(f'wbf_cover_>=2/{R_PER_IMAGE}', ascending=False)\nprint(f'WBF retention and coverage at IoU={WBF_IOU_THR}, R={R_PER_IMAGE}:')\nprint(wbf_analysis_df.to_string(index=False, float_format='%.3f'))\n\nfig, ax = plt.subplots(figsize=(12, 5))\nplot_df = wbf_analysis_df.sort_values(f'wbf_cover_>=2/{R_PER_IMAGE}')\nax.barh(plot_df['class_name'], plot_df[f'wbf_cover_>=2/{R_PER_IMAGE}'] * 100.0, color='#2a9d8f', alpha=0.8)\nax.set_xlabel(f'Image coverage with WBF score >= 2/{R_PER_IMAGE} (%)')\nax.set_title('WBF consensus coverage by class')\nplt.tight_layout()\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:44:53.779396Z","iopub.execute_input":"2026-06-01T13:44:53.779835Z","iopub.status.idle":"2026-06-01T13:44:54.146708Z","shell.execute_reply.started":"2026-06-01T13:44:53.779799Z","shell.execute_reply":"2026-06-01T13:44:54.145666Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 10. Co-occurrence Analysis\nMulti-label co-occurrence patterns reveal which abnormalities tend to appear together, informing classifier head design and data augmentation strategies.","metadata":{}},{"cell_type":"code","source":"# 10a. Co-occurrence heatmap\nimg_labels = train_df[train_df['class_id'] != 14].drop_duplicates(subset=['image_id', 'class_id'])\nimg_label_matrix = img_labels.pivot_table(index='image_id', columns='class_id', aggfunc='size', fill_value=0)\nimg_label_matrix = (img_label_matrix > 0).astype(int)\n\n# Co-occurrence matrix\ncooccur = img_label_matrix.T.dot(img_label_matrix)\nnp.fill_diagonal(cooccur.values, 0)\n\n# Normalize to conditional probability\nclass_freq = img_label_matrix.sum()\ncooccur_norm = cooccur.div(class_freq, axis=0)\n\nfig, ax = plt.subplots(figsize=(12, 10))\nlabels = [CLASS_MAP[i] for i in range(14)]\nsns.heatmap(cooccur_norm, ax=ax, xticklabels=labels, yticklabels=labels,\n            cmap='YlOrRd', vmin=0, fmt='.2f', annot=True, annot_kws={'size': 7})\nax.set_title('Conditional Co-occurrence P(class_j | class_i)', fontsize=13)\nplt.xticks(rotation=45, ha='right')\nplt.yticks(rotation=0)\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:44:54.14796Z","iopub.execute_input":"2026-06-01T13:44:54.148255Z","iopub.status.idle":"2026-06-01T13:44:54.941003Z","shell.execute_reply.started":"2026-06-01T13:44:54.148227Z","shell.execute_reply":"2026-06-01T13:44:54.939951Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 10b. Multi-label distribution + top pairs\nlabel_counts = img_label_matrix.sum(axis=1)\nprint('Number of abnormal labels per image:')\nprint(label_counts.value_counts().sort_index().to_string())\nprint(f'\\nMean labels per abnormal image: {label_counts.mean():.2f}')\nprint(f'Images with >= 3 labels: {(label_counts >= 3).sum()} ({(label_counts >= 3).mean():.1%})')\n\n# Top co-occurring pairs\npairs = []\nfor i, j in combinations(range(14), 2):\n    n_both = ((img_label_matrix[i] == 1) & (img_label_matrix[j] == 1)).sum()\n    if n_both > 0:\n        pairs.append({'class_i': CLASS_MAP[i], 'class_j': CLASS_MAP[j], 'n_cooccur': n_both})\n\npairs_df = pd.DataFrame(pairs).sort_values('n_cooccur', ascending=False)\nprint(f'\\nTop 10 co-occurring pairs:')\nprint(pairs_df.head(10).to_string(index=False))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:44:54.942229Z","iopub.execute_input":"2026-06-01T13:44:54.942579Z","iopub.status.idle":"2026-06-01T13:44:54.985148Z","shell.execute_reply.started":"2026-06-01T13:44:54.942521Z","shell.execute_reply":"2026-06-01T13:44:54.984134Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 11. Patient-Level Analysis\nCheck whether patient IDs are available for longitudinal analysis, which affects study design and train/val splitting strategy.","metadata":{}},{"cell_type":"code","source":"# 11. Patient ID availability\nif 'patient_id' in train_df.columns:\n    n_patients = train_df['patient_id'].nunique()\n    n_images = train_df['image_id'].nunique()\n    imgs_per_patient = train_df.groupby('patient_id')['image_id'].nunique()\n    print(f'Patient IDs available: Yes')\n    print(f'  Patients: {n_patients}, Images: {n_images}')\n    print(f'  Images per patient: mean={imgs_per_patient.mean():.2f}, max={imgs_per_patient.max()}')\n    print(f'  Patients with >1 image: {(imgs_per_patient > 1).sum()}')\n    print('  -> Use patient-level splitting to avoid data leakage.')\nelse:\n    print('Patient IDs not available in this dataset version.')\n    print('-> Cannot perform patient-level splitting.')\n    print('-> Must assume images are independent (potential leakage risk).')\n    print('   Recommendation: If patient IDs become available, re-split at patient level.')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:44:54.986392Z","iopub.execute_input":"2026-06-01T13:44:54.986784Z","iopub.status.idle":"2026-06-01T13:44:54.995115Z","shell.execute_reply.started":"2026-06-01T13:44:54.98675Z","shell.execute_reply":"2026-06-01T13:44:54.993314Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 12. Design Decisions Summary\nConsolidated table of all design recommendations derived from this EDA.","metadata":{}},{"cell_type":"code","source":"\n# 12. Design Decisions Summary\nwbf_cover_col = f'wbf_cover_>=2/{R_PER_IMAGE}'\ndesign_table = pd.DataFrame([\n    {'Component': 'Label Format', 'Decision': 'Soft labels (p_soft = n_positive / R_PER_IMAGE)', 'Rationale': f'{R_PER_IMAGE} radiologists/image from pool {R_POOL}; {(soft_summary_df[\"disagree_pct\"] > 0.3).sum()}/14 classes with >30% disagreement'},\n    {'Component': 'Loss Function', 'Decision': 'Quality Focal Loss with soft targets', 'Rationale': 'Class imbalance + annotation uncertainty'},\n    {'Component': 'IoU Threshold (theta_IoU)', 'Decision': 'Class-adaptive from symmetric best-match IoU', 'Rationale': f'Median IoU ranges from {iou_stats[\"median\"].min():.2f} to {iou_stats[\"median\"].max():.2f}'},\n    {'Component': 'WBF Consensus', 'Decision': 'Audit coverage at >=1/R, >=2/R, >=3/R', 'Rationale': f'Median >=2/R coverage: {wbf_di_df[wbf_cover_col].median():.2%}'},\n    {'Component': 'Input Resolution', 'Decision': '1024x1024 baseline for Stage 1/2; DINO exception 1022', 'Rationale': 'Matches current training notebooks and preserves more small-lesion detail than 512'},\n    {'Component': 'Feature Pyramid Strides', 'Decision': '[8, 16, 32, 64]', 'Rationale': 'Stride 8/16 remains important for focal classes and small boxes'},\n    {'Component': 'Detector Routing', 'Decision': 'Use spec DI formula: IoU_norm + (1-CV_area) + prevalence_norm', 'Rationale': 'DI is a routing policy from localization agreement, box geometry stability, and class support; disease probability remains Stage 1 p_cal'},\n    {'Component': 'Data Splitting', 'Decision': 'Image-level stratified unless patient IDs become available', 'Rationale': 'VinDr train.csv does not expose patient_id'},\n    {'Component': 'Preprocessing', 'Decision': 'VOI LUT + MONOCHROME1 fix', 'Rationale': 'DICOM metadata and training specs require robust CXR normalization'},\n    {'Component': 'Multi-label Head', 'Decision': 'Sigmoid (not Softmax)', 'Rationale': f'{(label_counts >= 2).sum()} images with >=2 concurrent abnormalities'},\n])\n\nprint(design_table.to_string(index=False))\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:44:54.996841Z","iopub.execute_input":"2026-06-01T13:44:54.997342Z","iopub.status.idle":"2026-06-01T13:44:55.024869Z","shell.execute_reply.started":"2026-06-01T13:44:54.997298Z","shell.execute_reply":"2026-06-01T13:44:55.023713Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"\n## 13. EDA DI Gate Audit\n\nDetectability Index (DI) is a class-level routing policy for the hybrid pipeline. It does **not** decide whether a disease is present; Stage 1 calibrated probability does that later. DI answers a narrower question: when Stage 1 is positive for a class, how suitable is that class for detector-based localization/refinement versus weak localization or deferral?\n\nThis cell follows the System Spec formula exactly:\n\n`DI(c) = w1 * IoU_norm(c) + w2 * (1 - CV_area(c)) + w3 * prevalence_norm(c)`\n\nThe three components are computed from earlier EDA tables only:\n\n- `IoU_norm(c)`: median symmetric best-match inter-radiologist IoU, already bounded to `[0, 1]`.\n- `CV_area(c)`: coefficient of variation of box area, using normalized `area_ratio` when image metadata is available.\n- `prevalence_norm(c)`: class positive-image count divided by the maximum positive-image count across the 14 disease classes.\n\nWeights default to equal `1/3` because the spec only constrains `w1 + w2 + w3 = 1`; any learned/tuned weights must be introduced later with a separate experiment. Route thresholds are operational audit thresholds for the baseline policy, not part of the DI formula and not a clinical metric.\n","metadata":{}},{"cell_type":"code","source":"# 13. EDA DI Gate Audit v3 - spec-faithful formula\n# DI is a routing heuristic, not disease probability and not a clinical metric.\n# System Spec formula:\n#   DI(c) = w1 * IoU_norm(c) + w2 * (1 - CV_area(c)) + w3 * prevalence_norm(c)\nimport json\nfrom pathlib import Path\n\nrequired_vars = ['train_df', 'bbox_df', 'iou_stats', 'CLASS_MAP']\nmissing_vars = [name for name in required_vars if name not in globals()]\nif missing_vars:\n    raise RuntimeError('DI gate audit requires running earlier EDA cells first. Missing: ' + ', '.join(missing_vars))\n\nif 'pd' not in globals() or 'np' not in globals():\n    raise RuntimeError('DI gate audit requires pandas as pd and numpy as np from the setup cell.')\n\nDI_SCORE_FORMULA = 'DI(c) = w1 * IoU_norm(c) + w2 * (1 - CV_area(c)) + w3 * prevalence_norm(c)'\nDI_WEIGHTS = {\n    'iou_norm': 1.0 / 3.0,\n    'one_minus_cv_area': 1.0 / 3.0,\n    'prevalence_norm': 1.0 / 3.0,\n}\nweight_sum = sum(float(value) for value in DI_WEIGHTS.values())\nif not np.isclose(weight_sum, 1.0):\n    raise RuntimeError(f'DI_WEIGHTS must sum to 1.0, got {weight_sum:.6f}')\n\n# These thresholds are policy/audit defaults, not part of the DI formula.\n# They must be tuned later with Stage 2 validation + Stage 3 risk analysis.\nDI_ROUTE_THRESHOLDS = {\n    'low_route_di': 0.35,\n    'high_route_di': 0.50,\n    'strong_positive_margin': 0.10,\n    'min_detector_consensus_boxes': 50,\n}\n\ndef _require_columns(df, cols, name):\n    missing = [col for col in cols if col not in df.columns]\n    if missing:\n        raise RuntimeError(f'{name} is missing required columns: {missing}')\n\ndef _to_float_or_none(value):\n    if value is None:\n        return None\n    try:\n        out = float(value)\n    except (TypeError, ValueError):\n        return None\n    if not np.isfinite(out):\n        return None\n    return out\n\ndef _round_or_none(value, ndigits=6):\n    out = _to_float_or_none(value)\n    return None if out is None else round(out, ndigits)\n\ndef _clip01_series(series):\n    values = pd.to_numeric(series, errors='coerce').replace([np.inf, -np.inf], np.nan)\n    return values.clip(0.0, 1.0)\n\ndef _route_tier(row):\n    if not bool(row['detector_supported']):\n        return 'not_detector_supported'\n    di_score = _to_float_or_none(row['di_score'])\n    if di_score is None:\n        return 'weak_preferred'\n    if di_score >= DI_ROUTE_THRESHOLDS['high_route_di']:\n        return 'detector_default'\n    if di_score >= DI_ROUTE_THRESHOLDS['low_route_di']:\n        return 'detector_if_strong_positive'\n    return 'weak_preferred'\n\ndef _resolve_project_output_root():\n    if Path('/kaggle').exists():\n        return Path('/kaggle/working/vindr_eda_di_gate')\n    cwd = Path.cwd().resolve()\n    if (cwd / 'artifacts').exists():\n        project_root = cwd\n    elif (cwd.parent / 'artifacts').exists():\n        project_root = cwd.parent\n    else:\n        project_root = cwd\n    return project_root / 'artifacts' / 'eda' / 'di_gate'\n\n# 1) IoU_norm(c): use the symmetric best-match median IoU from Section 8b.\niou_by_class = iou_stats.copy()\nif 'class_id' not in iou_by_class.columns:\n    iou_by_class = iou_by_class.reset_index().rename(columns={'index': 'class_id'})\n_require_columns(iou_by_class, ['class_id', 'median'], 'iou_stats')\niou_by_class = iou_by_class.set_index('class_id')\nmedian_iou = pd.to_numeric(iou_by_class['median'], errors='coerce').reindex(range(14)).fillna(0.0)\niou_norm = _clip01_series(median_iou).fillna(0.0)\n\n# 2) CV_area(c): coefficient of variation of box area.\n# Prefer normalized area_ratio so image-size differences do not dominate.\n_require_columns(bbox_df, ['class_id', 'image_id', 'area'], 'bbox_df')\narea_col = 'area_ratio' if 'area_ratio' in bbox_df.columns and bbox_df['area_ratio'].notna().any() else 'area'\narea_source = 'area_ratio' if area_col == 'area_ratio' else 'raw_pixel_area'\narea_df = bbox_df[['class_id', area_col]].copy()\narea_df[area_col] = pd.to_numeric(area_df[area_col], errors='coerce')\narea_df = area_df.replace([np.inf, -np.inf], np.nan).dropna(subset=[area_col])\narea_df = area_df[area_df[area_col] > 0]\nif area_df.empty:\n    raise RuntimeError('No positive bbox area values are available to compute CV_area(c).')\n\narea_stats = area_df.groupby('class_id')[area_col].agg(\n    n_area_values='count',\n    area_mean='mean',\n    area_std=lambda values: values.std(ddof=1) if len(values.dropna()) >= 2 else np.nan,\n).reindex(range(14))\ncv_area = area_stats['area_std'] / area_stats['area_mean']\ncv_area = cv_area.replace([np.inf, -np.inf], np.nan)\none_minus_cv_area = (1.0 - cv_area).clip(0.0, 1.0).fillna(0.0)\n\n# 3) prevalence_norm(c): positive-image count normalized across disease classes.\npositive_img_class = train_df[train_df['class_id'] != 14].drop_duplicates(['image_id', 'class_id'])\npositive_images = positive_img_class.groupby('class_id')['image_id'].nunique().reindex(range(14), fill_value=0).astype(float)\nmax_positive_images = float(positive_images.max())\nif max_positive_images <= 0:\n    raise RuntimeError('No positive disease images found; cannot compute prevalence_norm(c).')\nprevalence_norm = (positive_images / max_positive_images).clip(0.0, 1.0)\n\n# Detector support is an operational guard for routing. It is not a DI component.\n# Prefer WBF consensus boxes when the WBF audit table exists; otherwise fall back to raw bbox count.\nsupport_box_source = 'raw_boxes'\nsupport_boxes = bbox_df.groupby('class_id').size().reindex(range(14), fill_value=0).astype(int)\nif 'wbf_di_df' in globals():\n    wbf_tmp = wbf_di_df.copy()\n    if 'class_id' in wbf_tmp.columns:\n        wbf_tmp = wbf_tmp.set_index('class_id')\n        r_value = globals().get('R_PER_IMAGE', globals().get('R', None))\n        consensus_candidates = []\n        if r_value is not None:\n            consensus_candidates.extend([\n                f'wbf_boxes_>=2/{int(r_value)}',\n                f'wbf_boxes_>=1/{int(r_value)}',\n            ])\n        consensus_candidates.extend([col for col in wbf_tmp.columns if str(col).startswith('wbf_boxes_')])\n        consensus_col = next((col for col in consensus_candidates if col in wbf_tmp.columns), None)\n        if consensus_col is not None:\n            support_boxes = pd.to_numeric(wbf_tmp[consensus_col], errors='coerce').reindex(range(14), fill_value=0).fillna(0).astype(int)\n            support_box_source = consensus_col\n\ndetector_supported = support_boxes >= int(DI_ROUTE_THRESHOLDS['min_detector_consensus_boxes'])\n\nrecords = []\nfor cid in range(14):\n    class_name = CLASS_MAP[cid]\n    di_score = (\n        DI_WEIGHTS['iou_norm'] * float(iou_norm.loc[cid])\n        + DI_WEIGHTS['one_minus_cv_area'] * float(one_minus_cv_area.loc[cid])\n        + DI_WEIGHTS['prevalence_norm'] * float(prevalence_norm.loc[cid])\n    )\n    records.append({\n        'class_id': int(cid),\n        'class_name': class_name,\n        'iou_norm': float(iou_norm.loc[cid]),\n        'one_minus_cv_area': float(one_minus_cv_area.loc[cid]),\n        'prevalence_norm': float(prevalence_norm.loc[cid]),\n        'di_score': float(np.clip(di_score, 0.0, 1.0)),\n        'route_tier': None,\n        'detector_supported': bool(detector_supported.loc[cid]),\n        'n_positive_images': int(positive_images.loc[cid]),\n        'n_detector_support_boxes': int(support_boxes.loc[cid]),\n        'support_box_source': support_box_source,\n        'median_iou': float(median_iou.loc[cid]),\n        'cv_area': _to_float_or_none(cv_area.loc[cid]),\n        'area_mean': _to_float_or_none(area_stats.loc[cid, 'area_mean']) if cid in area_stats.index else None,\n        'area_std': _to_float_or_none(area_stats.loc[cid, 'area_std']) if cid in area_stats.index else None,\n        'n_area_values': int(area_stats.loc[cid, 'n_area_values']) if cid in area_stats.index and pd.notna(area_stats.loc[cid, 'n_area_values']) else 0,\n        'area_source': area_source,\n    })\n\ncomponents = pd.DataFrame(records)\ncomponents['route_tier'] = components.apply(_route_tier, axis=1)\ncomponents = components.sort_values('di_score', ascending=False).reset_index(drop=True)\n\npolicy_classes = {}\nfor row in components.sort_values('class_id').to_dict(orient='records'):\n    policy_classes[row['class_name']] = {\n        'class_id': int(row['class_id']),\n        'detector_supported': bool(row['detector_supported']),\n        'di_score': round(float(row['di_score']), 6),\n        'route_tier': row['route_tier'],\n        'components': {\n            'iou_norm': round(float(row['iou_norm']), 6),\n            'one_minus_cv_area': round(float(row['one_minus_cv_area']), 6),\n            'prevalence_norm': round(float(row['prevalence_norm']), 6),\n        },\n        'source_metrics': {\n            'median_iou': round(float(row['median_iou']), 6),\n            'cv_area': _round_or_none(row['cv_area']),\n            'area_mean': _round_or_none(row['area_mean']),\n            'area_std': _round_or_none(row['area_std']),\n            'area_source': row['area_source'],\n            'n_area_values': int(row['n_area_values']),\n            'n_positive_images': int(row['n_positive_images']),\n            'n_detector_support_boxes': int(row['n_detector_support_boxes']),\n            'support_box_source': row['support_box_source'],\n        },\n    }\n\npolicy = {\n    'policy_version': 'eda_di_gate_v3_spec_formula',\n    'policy_source': 'eda_di_audit_spec_v6_2_formula',\n    'policy_note': (\n        'DI uses the System Spec formula exactly: '\n        'DI(c)=w1*IoU_norm(c)+w2*(1-CV_area(c))+w3*prevalence_norm(c). '\n        'DI is a routing heuristic for Stage 1 positive records; it is not disease probability, '\n        'not detector mAP, and not a clinical metric.'\n    ),\n    'formula': DI_SCORE_FORMULA,\n    'normalization': {\n        'iou_norm': 'median symmetric best-match inter-radiologist IoU clipped to [0, 1]',\n        'cv_area': 'sample coefficient of variation std(area)/mean(area), using area_ratio when available',\n        'one_minus_cv_area': 'clip(1 - CV_area, 0, 1); undefined CV is conservatively scored as 0',\n        'prevalence_norm': 'n_positive_images / max_class_n_positive_images across the 14 disease classes',\n    },\n    'weights': DI_WEIGHTS,\n    'thresholds': DI_ROUTE_THRESHOLDS,\n    'threshold_note': (\n        'Route thresholds are operational audit defaults for the baseline and must be tuned with '\n        'Stage 2 validation and Stage 3 risk analysis. They are not part of the DI formula.'\n    ),\n    'route_logic': {\n        'detector_default': 'detector_supported and di_score >= high_route_di',\n        'detector_if_strong_positive': 'detector_supported and low_route_di <= di_score < high_route_di',\n        'weak_preferred': 'detector_supported but di_score < low_route_di, or DI undefined',\n        'not_detector_supported': 'support boxes below min_detector_consensus_boxes',\n    },\n    'classes': policy_classes,\n}\n\noutput_root = _resolve_project_output_root()\noutput_root.mkdir(parents=True, exist_ok=True)\ncomponents.to_csv(output_root / 'di_gate_audit.csv', index=False)\nwith (output_root / 'di_gate_policy.json').open('w', encoding='utf-8') as f:\n    json.dump(policy, f, indent=2, ensure_ascii=False, allow_nan=False)\n\nprint('DI gate audit exported:')\nprint(' -', output_root / 'di_gate_audit.csv')\nprint(' -', output_root / 'di_gate_policy.json')\nprint('\\nDI formula:', DI_SCORE_FORMULA)\nprint('DI weights:', DI_WEIGHTS)\nprint('Route thresholds are operational defaults, not formula terms:', DI_ROUTE_THRESHOLDS)\nprint('\\nDI routing table:')\nprint(components[[\n    'class_id', 'class_name', 'di_score', 'route_tier',\n    'iou_norm', 'one_minus_cv_area', 'prevalence_norm',\n    'cv_area', 'n_positive_images', 'n_detector_support_boxes'\n]].to_string(index=False, float_format='%.6f'))\n\nroute_counts = components['route_tier'].value_counts().rename_axis('route_tier').reset_index(name='n_classes')\nprint('\\nRoute distribution:')\nprint(route_counts.to_string(index=False))\n\nfig, ax = plt.subplots(figsize=(10, 5))\nplot_df = components.sort_values('di_score', ascending=True)\ncolors = plot_df['route_tier'].map({\n    'detector_default': '#2a9d8f',\n    'detector_if_strong_positive': '#e9c46a',\n    'weak_preferred': '#e76f51',\n    'not_detector_supported': '#8d99ae',\n}).fillna('#8d99ae')\nax.barh(plot_df['class_name'], plot_df['di_score'], color=colors)\nax.axvline(DI_ROUTE_THRESHOLDS['low_route_di'], color='#e9c46a', linestyle='--', linewidth=1, label='low_route_di')\nax.axvline(DI_ROUTE_THRESHOLDS['high_route_di'], color='#2a9d8f', linestyle='--', linewidth=1, label='high_route_di')\nax.set_xlabel('DI score')\nax.set_title('DI routing by class - spec formula')\nax.set_xlim(0, 1)\nax.legend(loc='lower right')\nplt.tight_layout()\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T13:44:55.026422Z","iopub.execute_input":"2026-06-01T13:44:55.027044Z","iopub.status.idle":"2026-06-01T13:44:55.396892Z","shell.execute_reply.started":"2026-06-01T13:44:55.026818Z","shell.execute_reply":"2026-06-01T13:44:55.39566Z"}},"outputs":[],"execution_count":null}]}