{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"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":[{"id":"2ae1e589-bd75-40a5-afce-70d1c259dc38","cell_type":"markdown","source":"# Interpretable Medical Anomaly Detection using DPM\n## Enhanced Notebook with Step-by-Step Visualizations & Comprehensive Metrics\n\n> **Dataset**: VinBigData Chest X-ray  \n> **Platform**: Kaggle Notebooks  \n> **Framework**: scikit-learn + OpenCV  \n\n### What this notebook adds over the base version:\n1. **Rich visualization at every pipeline step**\n2. **Anomaly summary output** — states exactly what was detected\n3. **DPM-paper-style annotations** — root filter (red) + part filters (colored dashed) + deformation ellipses\n4. **Root & part node decomposition** shown on X-ray images\n5. **Comprehensive metrics dashboard** — AP, Precision, Recall, FROC, F1/F2, Hinge Loss, Root/Part scores, Deformation costs\n","metadata":{}},{"id":"6c4367e1-d650-4786-837f-c7131ea32df1","cell_type":"markdown","source":"## Phase 0 — Environment Setup","metadata":{}},{"id":"79fbb960-d1b6-475a-877c-202a02e052d9","cell_type":"code","source":"!pip install -q pydicom opencv-python-headless scikit-learn scikit-image joblib tqdm","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"869ec81b-212e-4808-a15c-1eeef5fe70ef","cell_type":"code","source":"import os, sys, warnings, time\nimport numpy as np\nimport pandas as pd\nimport cv2\nimport pydicom\nimport matplotlib.pyplot as plt\nimport matplotlib.patches as mpatches\nfrom matplotlib.patches import Ellipse\nimport seaborn as sns\nfrom pathlib import Path\nfrom collections import defaultdict\nfrom tqdm import tqdm\nimport joblib\n\nfrom sklearn.svm import LinearSVC\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.metrics import classification_report\nfrom skimage.feature import hog\n\nwarnings.filterwarnings('ignore')\nnp.random.seed(42)\nsns.set_style('whitegrid')\nplt.rcParams.update({'figure.dpi': 120, 'font.size': 10})\nprint('Imports OK')","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"e3a1c2fb-36c7-48d8-b68d-cc77d3e04cb3","cell_type":"code","source":"from pathlib import Path\nBASE_DIR  = Path('/kaggle/input/competitions/vinbigdata-chest-xray-abnormalities-detection')\nOUTPUT_DIR = Path('/kaggle/working/output')\nOUTPUT_DIR.mkdir(parents=True, exist_ok=True)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"8407c80a-3cbc-4bf1-9b09-4f4420150fb6","cell_type":"code","source":"# ── Kaggle paths ──\nBASE_DIR  = Path('/kaggle/input/competitions/vinbigdata-chest-xray-abnormalities-detection')\nTRAIN_DIR = BASE_DIR / 'train'\nTEST_DIR  = BASE_DIR / 'test'\nTRAIN_CSV = BASE_DIR / 'train.csv'\nTEST_CSV  = BASE_DIR / 'test.csv'\n\nOUTPUT_DIR = Path('/kaggle/working/output')\nMODEL_DIR  = Path('/kaggle/working/models')\nOUTPUT_DIR.mkdir(parents=True, exist_ok=True)\nMODEL_DIR.mkdir(parents=True, exist_ok=True)\nprint('Paths set.')","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"4c7a7684-7df0-4a68-99d3-8d5508bf3d65","cell_type":"markdown","source":"## Phase 1 — Data Preparation & EDA\n\nVinDr-CXR has **multiple radiologist annotations per image**. We keep only findings agreed upon by ≥ 2 radiologists (consensus).\n","metadata":{}},{"id":"88bb11d9-3a00-4a02-902a-042f3e341b3c","cell_type":"code","source":"df = pd.read_csv(TRAIN_CSV)\nprint('Raw rows:', len(df))\nprint(df.head())","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"154610cc-f780-407b-abb6-842094b526e4","cell_type":"code","source":"# Separate findings from 'No finding'\ndf_findings = df[df['class_name'] != 'No finding'].copy()\ndf_normal   = df[df['class_name'] == 'No finding'].copy()\n\ndef merge_bboxes(group):\n    \"\"\"Union bounding box + count of radiologists who flagged it.\"\"\"\n    return pd.Series({\n        'x_min':  group['x_min'].min(),\n        'y_min':  group['y_min'].min(),\n        'x_max':  group['x_max'].max(),\n        'y_max':  group['y_max'].max(),\n        'count':  len(group)\n    })\n\ndf_merged    = df_findings.groupby(['image_id', 'class_name', 'class_id']).apply(merge_bboxes).reset_index()\ndf_consensus = df_merged[df_merged['count'] >= 2].copy()\nprint(f'Consensus annotations: {len(df_consensus)}')\nprint(df_consensus['class_name'].value_counts())","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"8fc2d631-af8f-4912-9600-63e47466df9c","cell_type":"code","source":"def load_dicom(image_id: str, split: str = 'train') -> np.ndarray:\n    \"\"\"Load a DICOM file and return a normalised uint8 grayscale array.\"\"\"\n    folder   = TRAIN_DIR if split == 'train' else TEST_DIR\n    dcm      = pydicom.dcmread(str(folder / f'{image_id}.dicom'))\n    img      = dcm.pixel_array.astype(np.float32)\n    # DICOM windowing\n    if hasattr(dcm, 'WindowCenter') and hasattr(dcm, 'WindowWidth'):\n        c = float(dcm.WindowCenter) if not isinstance(dcm.WindowCenter, pydicom.multival.MultiValue) else float(dcm.WindowCenter[0])\n        w = float(dcm.WindowWidth)  if not isinstance(dcm.WindowWidth,  pydicom.multival.MultiValue) else float(dcm.WindowWidth[0])\n        img = np.clip(img, c - w/2, c + w/2)\n    img = (img - img.min()) / (img.max() - img.min() + 1e-8) * 255\n    img = img.astype(np.uint8)\n    if dcm.PhotometricInterpretation == 'MONOCHROME1':\n        img = 255 - img\n    return img\n\n# Quick test\nsample_id = df_consensus['image_id'].iloc[3]\nsample_img = load_dicom(sample_id)\nprint(f'Sample image shape: {sample_img.shape}')","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"c7125f05-501f-4de4-80d4-47b5463da696","cell_type":"markdown","source":"### 📊 Visualization 1 — Exploratory Data Analysis","metadata":{}},{"id":"04b08dc7-4a41-4dc5-8d68-aedb8f26dd5a","cell_type":"code","source":"fig, axes = plt.subplots(2, 2, figsize=(16, 12))\n\n# (a) Class distribution\ncolors_bar = sns.color_palette('viridis', n_colors=len(df_consensus['class_name'].unique()))\ndf_consensus['class_name'].value_counts().plot(kind='barh', ax=axes[0,0], color=colors_bar)\naxes[0,0].set_title('Consensus Finding Distribution', fontweight='bold')\naxes[0,0].set_xlabel('Count')\n\n# (b) BBox size scatter colored by class\ndf_consensus['bbox_w'] = df_consensus['x_max'] - df_consensus['x_min']\ndf_consensus['bbox_h'] = df_consensus['y_max'] - df_consensus['y_min']\nscatter = axes[0,1].scatter(df_consensus['bbox_w'], df_consensus['bbox_h'],\n                            alpha=0.4, s=8, c=df_consensus['class_id'], cmap='tab10')\naxes[0,1].set_xlabel('BBox Width (px)'); axes[0,1].set_ylabel('BBox Height (px)')\naxes[0,1].set_title('Bounding Box Size by Class', fontweight='bold')\n\n# (c) Annotations per image histogram\nann_per_img = df_consensus.groupby('image_id').size()\naxes[1,0].hist(ann_per_img, bins=range(1, ann_per_img.max()+2), color='steelblue', edgecolor='white')\naxes[1,0].set_title('Annotations per Image', fontweight='bold')\naxes[1,0].set_xlabel('Number of annotations'); axes[1,0].set_ylabel('Frequency')\n\n# (d) BBox area distribution\ndf_consensus['bbox_area'] = df_consensus['bbox_w'] * df_consensus['bbox_h']\naxes[1,1].hist(df_consensus['bbox_area'], bins=50, color='coral', edgecolor='white')\naxes[1,1].set_title('Bounding Box Area Distribution', fontweight='bold')\naxes[1,1].set_xlabel('BBox Area (px²)')\n\nplt.suptitle('Phase 1 — EDA Overview', fontsize=15, fontweight='bold', y=1.01)\nplt.tight_layout()\nplt.savefig(OUTPUT_DIR / 'eda_overview.png', dpi=150, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"59192583-9d98-4da6-8847-945de8ab11d1","cell_type":"markdown","source":"### 📊 Visualization 2 — Sample GT Annotations (DPM-Paper-Style)","metadata":{}},{"id":"892ad805-72fd-43ba-9856-ede6e8770b87","cell_type":"code","source":"SHOW_CLASSES = ['Cardiomegaly', 'Pleural effusion', 'Pneumothorax', 'Aortic enlargement']\nfig, axes = plt.subplots(1, 4, figsize=(22, 5))\nfor ax, cls in zip(axes, SHOW_CLASSES):\n    subset = df_consensus[df_consensus['class_name'] == cls]\n    if subset.empty:\n        ax.set_visible(False); continue\n    row = subset.iloc[0]\n    img = load_dicom(row['image_id'])\n    ax.imshow(img, cmap='gray')\n    rect = plt.Rectangle((row['x_min'], row['y_min']),\n                          row['x_max'] - row['x_min'], row['y_max'] - row['y_min'],\n                          linewidth=3, edgecolor='red', facecolor='none')\n    ax.add_patch(rect)\n    cx = (row['x_min'] + row['x_max']) / 2\n    cy = (row['y_min'] + row['y_max']) / 2\n    ax.plot(cx, cy, 'r+', markersize=12, markeredgewidth=2)\n    ax.set_title(f'{cls}\\n({row[\"x_max\"]-row[\"x_min\"]:.0f}×{row[\"y_max\"]-row[\"y_min\"]:.0f} px)', fontsize=9)\n    ax.axis('off')\nplt.suptitle('Ground Truth — Root BBox Annotations (Red)', fontsize=13, fontweight='bold')\nplt.tight_layout()\nplt.savefig(OUTPUT_DIR / 'sample_annotations.png', dpi=150)\nplt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"84e49aa3-e3e8-4b8d-9fe9-14c43e0ba091","cell_type":"code","source":"TOP_CLASSES = ['Cardiomegaly', 'Pleural effusion', 'Aortic enlargement', 'Lung Opacity', 'Pneumothorax']\ndf_top      = df_consensus[df_consensus['class_name'].isin(TOP_CLASSES)].copy()\n\nunique_imgs = df_top['image_id'].unique()\nval_size    = int(0.2 * len(unique_imgs))\nval_ids     = set(np.random.choice(unique_imgs, val_size, replace=False))\ndf_top['split'] = df_top['image_id'].apply(lambda x: 'val' if x in val_ids else 'train')\nprint(df_top['split'].value_counts())","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"7929abe4-2d31-44e6-9ee0-f9f899a7ea0a","cell_type":"markdown","source":"## Phase 2 — Feature Extraction (HOG + Image Pyramid)\n","metadata":{}},{"id":"f71ca7d6-2cff-41b3-bb5e-d5ff7069f9b7","cell_type":"code","source":"class HOGExtractor:\n    \"\"\"HOG feature extractor with CLAHE preprocessing, matching DPM paper parameters.\"\"\"\n    def __init__(self, orientations=9, pixels_per_cell=(8, 8),\n                 cells_per_block=(2, 2), target_size=(512, 512)):\n        self.orientations    = orientations\n        self.pixels_per_cell = pixels_per_cell\n        self.cells_per_block = cells_per_block\n        self.target_size     = target_size\n\n    def preprocess(self, img: np.ndarray) -> np.ndarray:\n        img_r = cv2.resize(img, self.target_size)\n        clahe = cv2.createCLAHE(clipLimit=2.0, tileGridSize=(8, 8))\n        return clahe.apply(img_r)\n\n    def extract(self, img: np.ndarray):\n        \"\"\"Return (feature_vector, hog_image).\"\"\"\n        proc = self.preprocess(img)\n        return hog(proc, orientations=self.orientations,\n                   pixels_per_cell=self.pixels_per_cell,\n                   cells_per_block=self.cells_per_block,\n                   visualize=True, channel_axis=None)\n\nhog_ext = HOGExtractor()\nfeat_vec, hog_img = hog_ext.extract(sample_img)\nprint(f'HOG feature vector length: {len(feat_vec)}')","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"3c2d4983-c5d7-42b4-9e48-a036b850f985","cell_type":"markdown","source":"### 📊 Visualization 3 — HOG Feature Extraction","metadata":{}},{"id":"743038d4-0066-4817-ba58-f29f9b210cca","cell_type":"code","source":"fig, axes = plt.subplots(1, 3, figsize=(18, 5))\naxes[0].imshow(sample_img, cmap='gray')\naxes[0].set_title('Original CXR', fontweight='bold'); axes[0].axis('off')\n\naxes[1].imshow(hog_img, cmap='hot')\naxes[1].set_title('HOG Feature Map\\n(Edge Structure)', fontweight='bold'); axes[1].axis('off')\n\n# Overlay — shows HOG edges on anatomy\naxes[2].imshow(cv2.resize(sample_img, hog_img.shape[::-1]), cmap='gray', alpha=0.5)\naxes[2].imshow(hog_img, cmap='hot', alpha=0.5)\naxes[2].set_title('HOG Overlay on CXR', fontweight='bold'); axes[2].axis('off')\n\nplt.suptitle('Phase 2 — HOG Feature Extraction', fontsize=14, fontweight='bold')\nplt.tight_layout(); plt.savefig(OUTPUT_DIR / 'hog_visualisation.png', dpi=150); plt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"93806ee8-c305-4134-b5c6-c6caa757497a","cell_type":"code","source":"class ImagePyramid:\n    \"\"\"Gaussian image pyramid for multiscale DPM detection.\"\"\"\n    def __init__(self, n_scales=5, scale_factor=0.85, min_size=128):\n        self.n_scales     = n_scales\n        self.scale_factor = scale_factor\n        self.min_size     = min_size\n\n    def build(self, img: np.ndarray):\n        \"\"\"Returns list of (effective_scale, scaled_image) tuples.\"\"\"\n        pyramid = [(1.0, img)]\n        current = img.copy()\n        for lvl in range(1, self.n_scales):\n            h, w   = current.shape[:2]\n            new_h  = int(h * self.scale_factor)\n            new_w  = int(w * self.scale_factor)\n            if min(new_h, new_w) < self.min_size:\n                break\n            current = cv2.resize(current, (new_w, new_h))\n            pyramid.append((self.scale_factor ** lvl, current))\n        return pyramid","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"02483ece-fe10-4be3-a266-8e589fcf449f","cell_type":"markdown","source":"### 📊 Visualization 4 — Image Pyramid + HOG at Each Scale","metadata":{}},{"id":"e12cf4ff-0ff3-4a8f-b8e5-111d027fbc6d","cell_type":"code","source":"pyr = ImagePyramid().build(sample_img)\nfig, axes = plt.subplots(2, len(pyr), figsize=(18, 8))\nfor i, (scale, lvl_img) in enumerate(pyr):\n    axes[0, i].imshow(lvl_img, cmap='gray')\n    axes[0, i].set_title(f'Scale {scale:.2f}\\n{lvl_img.shape}', fontsize=8)\n    axes[0, i].axis('off')\n    _, hog_s = HOGExtractor(target_size=lvl_img.shape[::-1]).extract(lvl_img)\n    axes[1, i].imshow(hog_s, cmap='hot')\n    axes[1, i].set_title(f'HOG @ {scale:.2f}', fontsize=8)\n    axes[1, i].axis('off')\nfor j in range(len(pyr), axes.shape[1]):\n    axes[0, j].axis('off'); axes[1, j].axis('off')\nplt.suptitle('Image Pyramid with HOG at Each Scale', fontsize=13, fontweight='bold')\nplt.tight_layout(); plt.savefig(OUTPUT_DIR / 'image_pyramid.png', dpi=150); plt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"dbafdba6-4253-445b-9bd7-78124e4094e7","cell_type":"code","source":"def extract_patch_features(image_id, x_min, y_min, x_max, y_max, patch_size=(64, 64)):\n    \"\"\"Crop GT ROI and return HOG feature vector, or None if patch is empty.\"\"\"\n    img = load_dicom(image_id)\n    H, W = img.shape\n    x_min, y_min = max(0, int(x_min)), max(0, int(y_min))\n    x_max, y_max = min(W, int(x_max)), min(H, int(y_max))\n    patch = img[y_min:y_max, x_min:x_max]\n    if patch.size == 0:\n        return None\n    patch = cv2.resize(patch, patch_size)\n    feat, _ = hog(patch, orientations=9, pixels_per_cell=(8, 8),\n                  cells_per_block=(2, 2), visualize=True, channel_axis=None)\n    return feat\n\n\ndef extract_negative_patches(image_id, pos_boxes, patch_size=(64, 64), n_neg=5):\n    \"\"\"Sample random patches with IoU < 0.3 w.r.t. all GT boxes.\"\"\"\n    img = load_dicom(image_id)\n    H, W = img.shape\n    negatives, attempts = [], 0\n    while len(negatives) < n_neg and attempts < 100:\n        attempts += 1\n        px = np.random.randint(0, max(1, W - patch_size[0]))\n        py = np.random.randint(0, max(1, H - patch_size[1]))\n        iou_max = 0.0\n        for (x1, y1, x2, y2) in pos_boxes:\n            ix1, iy1 = max(px, x1), max(py, y1)\n            ix2, iy2 = min(px + patch_size[0], x2), min(py + patch_size[1], y2)\n            inter = max(0, ix2 - ix1) * max(0, iy2 - iy1)\n            union = patch_size[0] * patch_size[1] + (x2-x1)*(y2-y1) - inter\n            iou_max = max(iou_max, inter / (union + 1e-8))\n        if iou_max < 0.3:\n            p = img[py:py + patch_size[1], px:px + patch_size[0]]\n            f, _ = hog(p, orientations=9, pixels_per_cell=(8, 8),\n                       cells_per_block=(2, 2), visualize=True, channel_axis=None)\n            negatives.append(f)\n    return negatives","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"1f6b8866-28c7-4d20-8afe-adfaeb723593","cell_type":"markdown","source":"## Phase 3 — Root Filter Training (Global Detector)\n\nA **LinearSVC** is trained per class. Its weight vector `w` IS the root filter.\n","metadata":{}},{"id":"7080fba3-52a8-4c71-b779-e6388f868161","cell_type":"code","source":"TARGET_CLASS = 'Cardiomegaly'   # Change to loop over TOP_CLASSES later\ndf_train_cls = df_top[(df_top['class_name'] == TARGET_CLASS) & (df_top['split'] == 'train')]\nprint(f'Training positives for {TARGET_CLASS}: {len(df_train_cls)}')","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"9697ffe7-aab0-4e0d-b6b2-9031626b9f7a","cell_type":"code","source":"X_pos, X_neg = [], []\n\n# ── Positive patches ──\nfor _, row in tqdm(df_train_cls.iterrows(), total=len(df_train_cls), desc='Pos patches'):\n    feat = extract_patch_features(row['image_id'],\n                                   row['x_min'], row['y_min'],\n                                   row['x_max'], row['y_max'])\n    if feat is not None:\n        X_pos.append(feat)\n\n# ── Negative patches (from 'No finding' images) ──\nneg_img_ids = df_normal['image_id'].unique()[:100]\nfor img_id in tqdm(neg_img_ids, desc='Neg patches'):\n    negs = extract_negative_patches(img_id, [], n_neg=3)\n    X_neg.extend(negs)\n\nX = np.array(X_pos + X_neg)\ny = np.array([1]*len(X_pos) + [-1]*len(X_neg))\nprint(f'Positives: {len(X_pos)} | Negatives: {len(X_neg)} | Feature dim: {X.shape[1]}')","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"f27fe4a0-2f53-472a-8eb8-2939cc189cf3","cell_type":"code","source":"# ── Train root filter ──\nscaler   = StandardScaler()\nX_scaled = scaler.fit_transform(X)\n\nroot_svm = LinearSVC(C=0.01, max_iter=5000, class_weight='balanced')\nroot_svm.fit(X_scaled, y)\n\nw_root = root_svm.coef_[0]\nprint(f'Root filter shape: {w_root.shape}')\nprint(f'Root filter L2 norm: {np.linalg.norm(w_root):.4f}')\nprint(classification_report(y, root_svm.predict(X_scaled), target_names=['Negative', 'Positive']))","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"0785fbd9-2d91-4aed-b3d7-a2536ae7c014","cell_type":"code","source":"save_path = MODEL_DIR / f\"root_filter_{TARGET_CLASS.replace(' ', '_')}.pkl\"\njoblib.dump({'svm': root_svm, 'scaler': scaler, 'class': TARGET_CLASS}, save_path)\nprint(f'Saved: {save_path}')","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"7be13a6e-751b-49cb-bf75-8267c61b83c3","cell_type":"markdown","source":"### 📊 Visualization 5 — Root Filter Weights (DPM Paper Fig. 5 Style)","metadata":{}},{"id":"73ce86bf-00b6-49f1-8ee6-8a296d50ca6f","cell_type":"code","source":"fig, axes = plt.subplots(1, 3, figsize=(16, 5))\ntry:\n    w_img = w_root.reshape(7, 7, 36).mean(axis=2)\n    axes[0].imshow(w_img, cmap='RdBu_r', interpolation='nearest')\n    axes[0].set_title(f'Root Filter Weights\\n(Red=discriminative, Blue=suppressive)', fontweight='bold')\n    axes[0].axis('off')\nexcept Exception as e:\n    axes[0].hist(w_root, bins=50, color='steelblue')\n    axes[0].set_title('Root Filter Weight Distribution')\n\n# Weight histogram\naxes[1].hist(w_root, bins=50, color='steelblue', edgecolor='white', alpha=0.8)\naxes[1].axvline(0, color='red', linestyle='--', linewidth=1)\naxes[1].set_title('Root Filter Weight Distribution', fontweight='bold')\naxes[1].set_xlabel('Weight value')\n\n# Root filter score separation on training data\nscores_train = root_svm.decision_function(X_scaled)\naxes[2].hist(scores_train[y==1], bins=30, alpha=0.6, color='green', label='Positive', edgecolor='white')\naxes[2].hist(scores_train[y==-1], bins=30, alpha=0.6, color='red', label='Negative', edgecolor='white')\naxes[2].axvline(0, color='black', linestyle='--', linewidth=1)\naxes[2].set_title('Root Filter Score Distribution', fontweight='bold')\naxes[2].set_xlabel('SVM Decision Score'); axes[2].legend()\n\nplt.suptitle(f'Phase 3 — Root Filter Analysis: {TARGET_CLASS}', fontsize=14, fontweight='bold')\nplt.tight_layout(); plt.savefig(OUTPUT_DIR / f'root_filter_{TARGET_CLASS}.png', dpi=150); plt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"d7ac9228-cd0a-4116-869a-c830b437032e","cell_type":"markdown","source":"## Phase 4 — Part Filter Design & Latent SVM\n\nEach class has **anatomical part anchors** defined relative to the root bounding box.\n","metadata":{}},{"id":"e96afd6d-9baa-4732-bd86-f6f27630ccbf","cell_type":"code","source":"# ── Anatomical part anchors ──\nPART_ANCHORS = {\n    'Cardiomegaly': [\n        (0.0,  0.2,  0.3, 0.6),   # Left cardiac border\n        (0.7,  0.2,  0.3, 0.6),   # Right cardiac border\n        (0.1,  0.0,  0.8, 0.3),   # Superior contour\n        (0.1,  0.7,  0.8, 0.3),   # Inferior contour\n    ],\n    'Pleural effusion': [\n        (0.0,  0.6,  0.45, 0.4),\n        (0.55, 0.6,  0.45, 0.4),\n    ],\n    'Pneumothorax': [\n        (0.0,  0.0,  0.5,  0.4),\n        (0.5,  0.0,  0.5,  0.4),\n    ],\n    'Aortic enlargement': [\n        (0.2,  0.0,  0.6, 0.5),\n    ],\n    'Lung Opacity': [\n        (0.0,  0.0,  0.5, 0.5),\n        (0.5,  0.0,  0.5, 0.5),\n        (0.0,  0.5,  0.5, 0.5),\n        (0.5,  0.5,  0.5, 0.5),\n    ],\n}\nPART_NAMES = {\n    'Cardiomegaly': ['L-Border', 'R-Border', 'Superior', 'Inferior'],\n    'Pleural effusion': ['L-Costophrenic', 'R-Costophrenic'],\n    'Pneumothorax': ['L-Apex', 'R-Apex'],\n    'Aortic enlargement': ['Sup-Mediastinum'],\n    'Lung Opacity': ['UL-Quad', 'UR-Quad', 'LL-Quad', 'LR-Quad'],\n}\nPART_COLORS = ['cyan', 'lime', 'orange', 'magenta', 'yellow']\nprint('Part anchors defined for:', list(PART_ANCHORS.keys()))","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"b1718139-93be-4d3e-b76b-996682626927","cell_type":"code","source":"row = df_train_cls.iloc[0]\nimg = load_dicom(row['image_id'])\nanchors = PART_ANCHORS[TARGET_CLASS]\nnames   = PART_NAMES[TARGET_CLASS]\n\nfig, axes = plt.subplots(1, 3, figsize=(20, 6))\n\n# (a) Root only\naxes[0].imshow(img, cmap='gray')\nrect = plt.Rectangle((row['x_min'], row['y_min']),\n                      row['x_max']-row['x_min'], row['y_max']-row['y_min'],\n                      linewidth=3, edgecolor='red', facecolor='none')\naxes[0].add_patch(rect)\naxes[0].set_title('(a) Root Filter\\n(Global Detection)', fontweight='bold'); axes[0].axis('off')\n\n# (b) Parts only\naxes[1].imshow(img, cmap='gray')\nbw, bh = row['x_max']-row['x_min'], row['y_max']-row['y_min']\nfor i, (rx, ry, fw, fh) in enumerate(anchors):\n    px1 = int(row['x_min'] + rx * bw); py1 = int(row['y_min'] + ry * bh)\n    pw  = int(fw * bw); ph = int(fh * bh)\n    prect = plt.Rectangle((px1, py1), pw, ph,\n                           linewidth=2, edgecolor=PART_COLORS[i], facecolor='none', linestyle='--')\n    axes[1].add_patch(prect)\n    axes[1].text(px1+2, py1+12, f'P{i}: {names[i]}', color=PART_COLORS[i], fontsize=7,\n                 bbox=dict(boxstyle='round,pad=0.2', facecolor='black', alpha=0.6))\naxes[1].set_title('(b) Part Filters\\n(Anatomical Sub-Regions)', fontweight='bold'); axes[1].axis('off')\n\n# (c) Root + Parts combined (DPM paper style: red root, blue/colored parts)\naxes[2].imshow(img, cmap='gray')\nrect = plt.Rectangle((row['x_min'], row['y_min']),\n                      bw, bh, linewidth=3, edgecolor='red', facecolor='none')\naxes[2].add_patch(rect)\naxes[2].text(row['x_min'], row['y_min']-5, f'ROOT: {TARGET_CLASS}', color='red', fontsize=9,\n             bbox=dict(facecolor='white', alpha=0.8))\nfor i, (rx, ry, fw, fh) in enumerate(anchors):\n    px1 = int(row['x_min'] + rx * bw); py1 = int(row['y_min'] + ry * bh)\n    pw  = int(fw * bw); ph = int(fh * bh)\n    prect = plt.Rectangle((px1, py1), pw, ph, linewidth=2,\n                           edgecolor=PART_COLORS[i], facecolor=PART_COLORS[i],\n                           linestyle='--', alpha=0.15)\n    axes[2].add_patch(prect)\n    prect2 = plt.Rectangle((px1, py1), pw, ph, linewidth=2,\n                            edgecolor=PART_COLORS[i], facecolor='none', linestyle='--')\n    axes[2].add_patch(prect2)\n    pcx, pcy = px1 + pw//2, py1 + ph//2\n    axes[2].plot(pcx, pcy, 'o', color=PART_COLORS[i], markersize=6)\n    axes[2].text(px1+2, py1+12, f'P{i}', color=PART_COLORS[i], fontsize=8, fontweight='bold')\naxes[2].set_title('(c) Root + Parts Combined\\n(DPM Paper Style)', fontweight='bold'); axes[2].axis('off')\n\nplt.suptitle(f'DPM Model Structure — {TARGET_CLASS}\\nRed=Root, Colored=Parts', fontsize=14, fontweight='bold')\nplt.tight_layout(); plt.savefig(OUTPUT_DIR / 'root_parts_decomposition.png', dpi=150); plt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"b62561e5-9f64-4794-92d2-d0af7295646f","cell_type":"markdown","source":"### 📊 Visualization 6 — Root + Part Decomposition on X-ray (DPM Paper Fig. 2/7 Style)","metadata":{}},{"id":"7b8af668-52df-4a76-affe-0c50c9ed71e4","cell_type":"code","source":"def extract_part_patches(image_id, x_min, y_min, x_max, y_max,\n                          anchors, part_size=(32, 32)):\n    \"\"\"Extract HOG features for each anatomical part sub-region.\"\"\"\n    img = load_dicom(image_id)\n    H, W = img.shape\n    bw, bh = x_max - x_min, y_max - y_min\n    part_features = []\n    for (rx, ry, fw, fh) in anchors:\n        px1 = max(0, int(x_min + rx * bw))\n        py1 = max(0, int(y_min + ry * bh))\n        px2 = min(W, int(px1 + fw * bw))\n        py2 = min(H, int(py1 + fh * bh))\n        patch = img[py1:py2, px1:px2]\n        if patch.size == 0:\n            part_features.append(None)\n            continue\n        patch = cv2.resize(patch, part_size)\n        f, _  = hog(patch, orientations=9, pixels_per_cell=(4, 4),\n                    cells_per_block=(2, 2), visualize=True, channel_axis=None)\n        part_features.append(f)\n    return part_features","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"c43a2bae-1065-4221-8281-44d043a89238","cell_type":"code","source":"class LatentSVM:\n    \"\"\"\n    Simplified Latent SVM via coordinate descent.\n    Combines root features + part features into one concatenated feature vector.\n    \"\"\"\n    def __init__(self, n_parts, C=0.01, max_iter=5):\n        self.n_parts      = n_parts\n        self.C            = C\n        self.max_iter     = max_iter\n        self.combined_svm = None\n        self.scaler       = None\n        self.deform_w     = [np.array([0.01, 0.001, 0.01, 0.001])] * n_parts\n\n    def fit(self, X_root, X_parts_list, y):\n        valid_parts = [p for p in X_parts_list if p is not None]\n        X_all = np.hstack([X_root] + valid_parts)\n        self.scaler = StandardScaler()\n        X_scaled    = self.scaler.fit_transform(X_all)\n\n        for it in range(self.max_iter):\n            svm = LinearSVC(C=self.C, max_iter=3000, class_weight='balanced')\n            svm.fit(X_scaled, y)\n            self.combined_svm = svm\n            scores = svm.decision_function(X_scaled)\n            hard_neg = np.sum((y == -1) & (scores > -1.0))\n            print(f'  Iter {it+1}/{self.max_iter}  |  '\n                  f'Train acc: {svm.score(X_scaled, y):.4f}  |  '\n                  f'Hard negs: {hard_neg}')\n        return self\n\n    def decision_function(self, X_root, X_parts_list):\n        valid_parts = [p for p in X_parts_list if p is not None]\n        X_all    = np.hstack([X_root] + valid_parts)\n        X_scaled = self.scaler.transform(X_all)\n        return self.combined_svm.decision_function(X_scaled)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"d7a95cf5-b548-4184-8274-1306d59dd658","cell_type":"code","source":"# ── Build root + part feature matrices for Latent SVM ──\nanchors = PART_ANCHORS[TARGET_CLASS]\nn_parts = len(anchors)\n\nX_root_train, X_parts_train, y_lsvm = [], [[] for _ in range(n_parts)], []\n\nprint(f'Building LSVM features for {TARGET_CLASS} ({n_parts} parts)...')\nfor _, row in tqdm(df_train_cls.iterrows(), total=len(df_train_cls)):\n    rf = extract_patch_features(row['image_id'],\n                                 row['x_min'], row['y_min'],\n                                 row['x_max'], row['y_max'])\n    pf = extract_part_patches(row['image_id'],\n                               row['x_min'], row['y_min'],\n                               row['x_max'], row['y_max'], anchors)\n    if rf is None or any(p is None for p in pf):\n        continue\n    X_root_train.append(rf)\n    for i, p in enumerate(pf):\n        X_parts_train[i].append(p)\n    y_lsvm.append(1)\n\n# Add negatives (root only; parts filled with zeros)\ndummy_part = np.zeros_like(X_parts_train[0][0]) if X_parts_train[0] else None\nfor img_id in tqdm(neg_img_ids[:50], desc='Neg (LSVM)'):\n    negs = extract_negative_patches(img_id, [], n_neg=2)\n    neg_parts = [np.zeros_like(X_parts_train[i][0]) for i in range(n_parts)] \\\n                if X_parts_train[0] else []\n    for nf in negs:\n        X_root_train.append(nf)\n        for i in range(n_parts):\n            X_parts_train[i].append(neg_parts[i])\n        y_lsvm.append(-1)\n\nX_root_arr  = np.array(X_root_train)\nX_parts_arr = [np.array(X_parts_train[i]) for i in range(n_parts)]\ny_lsvm      = np.array(y_lsvm)\nprint(f'LSVM dataset: {X_root_arr.shape[0]} samples, root_dim={X_root_arr.shape[1]}')","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"2cd59db2-be5b-408b-af43-220bf30b5f43","cell_type":"code","source":"lsvm = LatentSVM(n_parts=n_parts, C=0.01, max_iter=5)\nlsvm.fit(X_root_arr, X_parts_arr, y_lsvm)\njoblib.dump(lsvm, MODEL_DIR / f'lsvm_{TARGET_CLASS.replace(\" \",\"_\")}.pkl')\nprint('Latent SVM saved.')","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"15b1494e-9ec3-44cd-9631-523d5b966958","cell_type":"markdown","source":"## Phase 5 — Deformation Model\n\nThe deformation cost uses the **quadratic spring model**: `cost_i = d_i · [dx, dx², dy, dy²]`\n","metadata":{}},{"id":"82e4a578-b976-42f0-ba89-86c201207935","cell_type":"code","source":"def deformation_score(actual_pos, anchor_pos, deform_weights):\n    \"\"\"Compute deformation penalty for one part. Returns a negative number (penalty).\"\"\"\n    dx    = actual_pos[0] - anchor_pos[0]\n    dy    = actual_pos[1] - anchor_pos[1]\n    phi_d = np.array([dx, dx**2, dy, dy**2], dtype=np.float32)\n    return -float(np.dot(deform_weights, phi_d))\n\ndef compute_full_dpm_score(root_score, part_scores, deform_scores):\n    \"\"\"Score = root_score + Σ(part_score_i + deformation_score_i)\"\"\"\n    return root_score + sum(ps + ds for ps, ds in zip(part_scores, deform_scores))\n\n# ── Quick sanity check ──\npenalty_on_anchor = deformation_score((100, 100), (100, 100), np.array([0.01, 0.001, 0.01, 0.001]))\npenalty_displaced = deformation_score((150, 130), (100, 100), np.array([0.01, 0.001, 0.01, 0.001]))\nprint(f'Deformation cost at anchor:    {penalty_on_anchor:.4f}  (should be 0)')\nprint(f'Deformation cost at +50,+30px: {penalty_displaced:.4f}  (negative penalty)')","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"d1193811-f6f4-4a3e-b762-0d59006609f0","cell_type":"markdown","source":"### 📊 Visualization 7 — Deformation Cost Surface","metadata":{}},{"id":"2d09d35f-c8fb-43af-8476-71bec3485f30","cell_type":"code","source":"fig, axes = plt.subplots(1, 2, figsize=(14, 5))\ndw = np.array([0.01, 0.001, 0.01, 0.001])\noffsets = np.arange(-60, 61, 2)\nX_grid, Y_grid = np.meshgrid(offsets, offsets)\nZ = np.zeros_like(X_grid, dtype=float)\nfor i in range(len(offsets)):\n    for j in range(len(offsets)):\n        Z[i, j] = deformation_score((100+offsets[j], 100+offsets[i]), (100, 100), dw)\n\nim = axes[0].imshow(Z, extent=[-60,60,-60,60], cmap='RdYlGn', origin='lower')\naxes[0].set_title('Deformation Cost Surface\\n(Green=low cost, Red=high penalty)', fontweight='bold')\naxes[0].set_xlabel('dx (pixels)'); axes[0].set_ylabel('dy (pixels)')\naxes[0].plot(0, 0, 'k+', markersize=15, markeredgewidth=2)\nplt.colorbar(im, ax=axes[0], label='Cost')\n\n# 1D slice\ncosts_1d = [deformation_score((100+d, 100), (100, 100), dw) for d in offsets]\naxes[1].plot(offsets, costs_1d, 'b-', linewidth=2)\naxes[1].axhline(0, color='gray', linestyle='--')\naxes[1].axvline(0, color='gray', linestyle='--')\naxes[1].set_title('Deformation Cost (1D slice, dy=0)', fontweight='bold')\naxes[1].set_xlabel('dx (pixels)'); axes[1].set_ylabel('Penalty')\naxes[1].fill_between(offsets, costs_1d, 0, alpha=0.2, color='blue')\n\nplt.suptitle('Phase 5 — Deformation Model (Quadratic Spring)', fontsize=14, fontweight='bold')\nplt.tight_layout(); plt.savefig(OUTPUT_DIR / 'deformation_cost.png', dpi=150); plt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"e2b4fcaf-d25c-41aa-bf08-1bb4e842ef53","cell_type":"markdown","source":"## Phase 6 — Sliding Window Inference + NMS\n","metadata":{}},{"id":"61d5bb59-474e-4106-8397-c2e7171a5e57","cell_type":"code","source":"def sliding_window_detect(img, svm_model, feat_scaler,\n                           patch_size=(64, 64), stride=16, threshold=0.3):\n    \"\"\"Multi-scale sliding window detector. Returns list of (score, x1, y1, x2, y2).\"\"\"\n    pyramid    = ImagePyramid(n_scales=5, scale_factor=0.85).build(img)\n    detections = []\n    for scale, scaled_img in pyramid:\n        H, W = scaled_img.shape\n        for row_y in range(0, H - patch_size[1], stride):\n            for col_x in range(0, W - patch_size[0], stride):\n                patch = scaled_img[row_y:row_y + patch_size[1],\n                                   col_x:col_x + patch_size[0]]\n                feat, _ = hog(patch, orientations=9, pixels_per_cell=(8, 8),\n                              cells_per_block=(2, 2), visualize=True, channel_axis=None)\n                feat_s  = feat_scaler.transform([feat])\n                score   = svm_model.decision_function(feat_s)[0]\n                if score > threshold:\n                    x1 = int(col_x / scale)\n                    y1 = int(row_y / scale)\n                    x2 = int((col_x + patch_size[0]) / scale)\n                    y2 = int((row_y + patch_size[1]) / scale)\n                    detections.append((score, x1, y1, x2, y2))\n    return detections\n\n\ndef nms(detections, iou_threshold=0.3):\n    \"\"\"Greedy Non-Maximum Suppression.\"\"\"\n    if not detections:\n        return []\n    dets = sorted(detections, key=lambda x: -x[0])\n    kept = []\n    while dets:\n        best = dets.pop(0)\n        kept.append(best)\n        _, bx1, by1, bx2, by2 = best\n        remaining = []\n        for d in dets:\n            _, dx1, dy1, dx2, dy2 = d\n            ix1, iy1 = max(bx1, dx1), max(by1, dy1)\n            ix2, iy2 = min(bx2, dx2), min(by2, dy2)\n            inter = max(0, ix2 - ix1) * max(0, iy2 - iy1)\n            union = (bx2-bx1)*(by2-by1) + (dx2-dx1)*(dy2-dy1) - inter\n            if inter / (union + 1e-8) < iou_threshold:\n                remaining.append(d)\n        dets = remaining\n    return kept","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"7564c598-a504-40fa-a8ca-040cc0c58695","cell_type":"code","source":"# ── Run detector on a validation image ──\nsaved_model = joblib.load(MODEL_DIR / f\"root_filter_{TARGET_CLASS.replace(' ', '_')}.pkl\")\nval_row  = df_top[(df_top['class_name'] == TARGET_CLASS) & (df_top['split'] == 'val')].iloc[0]\nval_img  = load_dicom(val_row['image_id'])\n\ndets = sliding_window_detect(val_img, saved_model['svm'], saved_model['scaler'],\n                              stride=32, threshold=0.2)\ndets = nms(dets, iou_threshold=0.3)\nprint(f'Detections after NMS: {len(dets)}')\nif dets:\n    print('Top-3 detections (score, x1, y1, x2, y2):')\n    for d in dets[:3]:\n        print(' ', d)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"46eed9c5-5e86-4ef4-b912-b05e9ea7329e","cell_type":"markdown","source":"### 📊 Visualization 8 — Sliding Window Detection + Before/After NMS","metadata":{}},{"id":"93707763-62de-44d9-93e5-2eabfdbefc50","cell_type":"code","source":"raw_dets = sliding_window_detect(val_img, saved_model['svm'], saved_model['scaler'], stride=32, threshold=0.1)\nnms_dets = nms(raw_dets, iou_threshold=0.3)\n\nfig, axes = plt.subplots(1, 3, figsize=(20, 6))\n# Before NMS\naxes[0].imshow(val_img, cmap='gray')\nfor s, x1, y1, x2, y2 in raw_dets[:50]:\n    rect = plt.Rectangle((x1, y1), x2-x1, y2-y1, linewidth=1, edgecolor='yellow', facecolor='none', alpha=0.4)\n    axes[0].add_patch(rect)\naxes[0].set_title(f'Before NMS\\n({len(raw_dets)} detections)', fontweight='bold'); axes[0].axis('off')\n\n# After NMS\naxes[1].imshow(val_img, cmap='gray')\nfor s, x1, y1, x2, y2 in nms_dets:\n    rect = plt.Rectangle((x1, y1), x2-x1, y2-y1, linewidth=3, edgecolor='red', facecolor='none')\n    axes[1].add_patch(rect)\n    axes[1].text(x1, y1-3, f'{s:.2f}', color='red', fontsize=8, bbox=dict(facecolor='white', alpha=0.7))\n# GT box (green dashed)\ngt_rect = plt.Rectangle((val_row['x_min'], val_row['y_min']),\n                         val_row['x_max']-val_row['x_min'], val_row['y_max']-val_row['y_min'],\n                         linewidth=2, edgecolor='green', facecolor='none', linestyle='--')\naxes[1].add_patch(gt_rect)\naxes[1].set_title(f'After NMS\\n({len(nms_dets)} kept, green=GT)', fontweight='bold'); axes[1].axis('off')\n\n# Root + parts on best detection\naxes[2].imshow(val_img, cmap='gray')\nif nms_dets:\n    s, x1, y1, x2, y2 = nms_dets[0]\n    rect = plt.Rectangle((x1, y1), x2-x1, y2-y1, linewidth=3, edgecolor='red', facecolor='none')\n    axes[2].add_patch(rect)\n    axes[2].text(x1, y1-5, f'ROOT: {s:.2f}', color='red', fontsize=9,\n                 bbox=dict(facecolor='white', alpha=0.8))\n    bw, bh = x2-x1, y2-y1\n    for i, (rx, ry, fw, fh) in enumerate(PART_ANCHORS[TARGET_CLASS]):\n        px1 = int(x1 + rx*bw); py1 = int(y1 + ry*bh)\n        pw = int(fw*bw); ph = int(fh*bh)\n        prect = plt.Rectangle((px1, py1), pw, ph, linewidth=2,\n                               edgecolor=PART_COLORS[i], facecolor='none', linestyle='--')\n        axes[2].add_patch(prect)\n        axes[2].text(px1+2, py1+10, f'P{i}', color=PART_COLORS[i], fontsize=8, fontweight='bold')\naxes[2].set_title('DPM Detection: Root + Parts', fontweight='bold'); axes[2].axis('off')\n\nplt.suptitle(f'Phase 6 — Sliding Window Detection ({TARGET_CLASS})', fontsize=14, fontweight='bold')\nplt.tight_layout(); plt.savefig(OUTPUT_DIR / 'detection_nms.png', dpi=150); plt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"bedab68d-5f3f-498d-8ce3-ec1bfc59ca8a","cell_type":"markdown","source":"## Phase 7 — Comprehensive Evaluation & Metrics Dashboard\n\nComputing: **AP, Precision, Recall, IoU, Hypothesis Score, Deformation Cost, Hinge Loss, FROC, F1/F2, Root Filter Score, Part Contribution Score**\n","metadata":{}},{"id":"ff6cc1de-1fea-41ff-ae78-c0c3fccc6984","cell_type":"code","source":"def compute_iou(b1, b2):\n    x1, y1 = max(b1[0], b2[0]), max(b1[1], b2[1])\n    x2, y2 = min(b1[2], b2[2]), min(b1[3], b2[3])\n    inter   = max(0, x2-x1) * max(0, y2-y1)\n    union   = (b1[2]-b1[0])*(b1[3]-b1[1]) + (b2[2]-b2[0])*(b2[3]-b2[1]) - inter\n    return inter / (union + 1e-8)\n\ndef compute_ap(preds, gt_dict, iou_thresh=0.5):\n    \"\"\"Returns (AP, recall_curve, precision_curve, tp_list, fp_list).\"\"\"\n    preds      = sorted(preds, key=lambda x: -x[0])\n    matched    = defaultdict(set)\n    tp, fp     = [], []\n    for score, x1, y1, x2, y2, img_id in preds:\n        gt_boxes = gt_dict.get(img_id, [])\n        hit      = False\n        for i, gb in enumerate(gt_boxes):\n            if i in matched[img_id]:\n                continue\n            if compute_iou((x1, y1, x2, y2), gb) >= iou_thresh:\n                tp.append(1); fp.append(0)\n                matched[img_id].add(i)\n                hit = True; break\n        if not hit:\n            tp.append(0); fp.append(1)\n    n_gt      = sum(len(v) for v in gt_dict.values())\n    tp_c      = np.cumsum(tp)\n    fp_c      = np.cumsum(fp)\n    recall    = tp_c / (n_gt + 1e-8)\n    precision = tp_c / (tp_c + fp_c + 1e-8)\n    return float(np.trapz(precision, recall)), recall, precision, tp, fp","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"d0d19421-fde2-4099-8cb3-a5fb6bb09a05","cell_type":"code","source":"# ── Evaluate on validation set ──\ndf_val_cls = df_top[(df_top['class_name'] == TARGET_CLASS) & (df_top['split'] == 'val')]\ngt_dict    = defaultdict(list)\nfor _, r in df_val_cls.iterrows():\n    gt_dict[r['image_id']].append((r['x_min'], r['y_min'], r['x_max'], r['y_max']))\n\nall_preds = []\nall_scores_list = []\nall_ious = []\nn_images_evaluated = 0\n\nfor img_id in tqdm(list(gt_dict.keys())[:30], desc='Evaluating'):\n    img  = load_dicom(img_id)\n    det = sliding_window_detect(img, saved_model['svm'], saved_model['scaler'],\n                                  stride=32, threshold=0.1)\n    det = nms(det)\n    n_images_evaluated += 1\n    for score, x1, y1, x2, y2 in det:\n        all_preds.append((score, x1, y1, x2, y2, img_id))\n        all_scores_list.append(score)\n        best_iou = max([compute_iou((x1,y1,x2,y2), gb) for gb in gt_dict[img_id]] or [0])\n        all_ious.append(best_iou)\n\nap_05, recall_curve, prec_curve, tp_list, fp_list = compute_ap(all_preds, gt_dict, iou_thresh=0.5)\nap_03, _, _, _, _ = compute_ap(all_preds, gt_dict, iou_thresh=0.3)\nprint(f'AP@0.5: {ap_05:.4f}')\nprint(f'AP@0.3: {ap_03:.4f}')","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"8c60fed3-3bcd-4318-b865-38f36c9ff7f9","cell_type":"code","source":" #── Compute ALL requested metrics ──\ntp_total = sum(tp_list)\nfp_total = sum(fp_list)\nn_gt_total = sum(len(v) for v in gt_dict.values())\nfn_total = n_gt_total - tp_total\n\nprecision_val = tp_total / (tp_total + fp_total + 1e-8)\nrecall_val    = tp_total / (n_gt_total + 1e-8)\nf1_val = 2 * precision_val * recall_val / (precision_val + recall_val + 1e-8)\nf2_val = 5 * precision_val * recall_val / (4 * precision_val + recall_val + 1e-8)\n\n# Hinge loss\ny_val_binary = np.array([1 if iou >= 0.5 else -1 for iou in all_ious])\nscores_arr = np.array(all_scores_list)\nhinge_losses = np.maximum(0, 1 - y_val_binary * scores_arr)\navg_hinge = np.mean(hinge_losses) if len(hinge_losses) > 0 else 0\n\n# Root filter scores\nroot_scores_pos = scores_arr[y_val_binary == 1] if np.any(y_val_binary == 1) else np.array([0])\nroot_scores_neg = scores_arr[y_val_binary == -1] if np.any(y_val_binary == -1) else np.array([0])\n\n# Part contribution\nif lsvm.combined_svm is not None:\n    w_combined = lsvm.combined_svm.coef_[0]\n    root_dim = X_root_arr.shape[1]\n    root_w_norm = np.linalg.norm(w_combined[:root_dim])\n    part_w_norm = np.linalg.norm(w_combined[root_dim:])\n    total_w_norm = root_w_norm + part_w_norm\n    root_contrib_pct = root_w_norm / total_w_norm * 100\n    part_contrib_pct = part_w_norm / total_w_norm * 100\nelse:\n    root_contrib_pct, part_contrib_pct = 50, 50\n\n# Deformation cost stats\ndeform_costs = []\nfor i in range(n_parts):\n    for d in range(10):\n        dc = deformation_score((100 + d*5, 100 + d*3), (100, 100), lsvm.deform_w[i])\n        deform_costs.append(dc)\navg_deform = np.mean(deform_costs)\navg_iou = np.mean(all_ious) if all_ious else 0\n\n# FROC curve\nfp_per_img_rates = [0.5, 1.0, 2.0, 4.0, 8.0]\nfroc_sensitivities = []\nsorted_preds = sorted(all_preds, key=lambda x: -x[0])\nfor fpi_target in fp_per_img_rates:\n    max_fp = int(fpi_target * n_images_evaluated)\n    matched_froc = defaultdict(set)\n    tp_f, fp_f = 0, 0\n    for score, x1, y1, x2, y2, img_id in sorted_preds:\n        if fp_f >= max_fp: break\n        hit = False\n        for i, gb in enumerate(gt_dict.get(img_id, [])):\n            if i in matched_froc[img_id]: continue\n            if compute_iou((x1,y1,x2,y2), gb) >= 0.5:\n                tp_f += 1; matched_froc[img_id].add(i); hit = True; break\n        if not hit:\n            fp_f += 1\n    froc_sensitivities.append(tp_f / (n_gt_total + 1e-8))\nprint('All metrics computed.')","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"8f396c1d-0b5b-4c11-88a5-218be9b37dc1","cell_type":"markdown","source":"### 📊 Visualization 9 — Metrics Dashboard","metadata":{}},{"id":"52fc6456-1de9-4778-8d4d-ee4a72279b03","cell_type":"code","source":"fig = plt.figure(figsize=(22, 18))\ngs = fig.add_gridspec(3, 3, hspace=0.35, wspace=0.3)\n\n# (a) Precision-Recall Curve\nax = fig.add_subplot(gs[0, 0])\nif len(recall_curve) > 0:\n    ax.plot(recall_curve, prec_curve, 'b-', linewidth=2)\n    ax.fill_between(recall_curve, prec_curve, alpha=0.2, color='blue')\nax.set_xlabel('Recall'); ax.set_ylabel('Precision')\nax.set_title(f'Precision-Recall Curve\\nAP@0.5 = {ap_05:.4f}', fontweight='bold')\nax.set_xlim([0, 1]); ax.set_ylim([0, 1])\n\n# (b) FROC Curve\nax = fig.add_subplot(gs[0, 1])\nax.plot(fp_per_img_rates, froc_sensitivities, 'ro-', linewidth=2, markersize=8)\nax.set_xlabel('False Positives per Image'); ax.set_ylabel('Sensitivity')\nax.set_title('FROC Curve', fontweight='bold')\nax.set_ylim([0, 1]); ax.grid(True, alpha=0.3)\nfor x, yv in zip(fp_per_img_rates, froc_sensitivities):\n    ax.annotate(f'{yv:.2f}', (x, yv), textcoords='offset points', xytext=(5, 10), fontsize=8)\n\n# (c) IoU Distribution (Bounding Box Overlap Threshold)\nax = fig.add_subplot(gs[0, 2])\nif all_ious:\n    ax.hist(all_ious, bins=20, color='teal', edgecolor='white', alpha=0.8)\n    ax.axvline(0.5, color='red', linestyle='--', label='IoU=0.5')\n    ax.axvline(0.3, color='orange', linestyle='--', label='IoU=0.3')\n    ax.legend(fontsize=8)\nax.set_xlabel('IoU'); ax.set_title(f'BBox Overlap Distribution\\nMean IoU={avg_iou:.3f}', fontweight='bold')\n\n# (d) F1/F2 + Precision/Recall bar chart\nax = fig.add_subplot(gs[1, 0])\nmetrics_names = ['Precision', 'Recall', 'F1 Score', 'F2 Score']\nmetrics_vals  = [precision_val, recall_val, f1_val, f2_val]\nbars = ax.barh(metrics_names, metrics_vals, color=['#3498db', '#2ecc71', '#e74c3c', '#9b59b6'])\nfor bar, val in zip(bars, metrics_vals):\n    ax.text(bar.get_width() + 0.02, bar.get_y() + bar.get_height()/2, f'{val:.3f}', va='center')\nax.set_xlim([0, 1.2]); ax.set_title('Classification Metrics', fontweight='bold')\n\n# (e) Hinge Loss\nax = fig.add_subplot(gs[1, 1])\nif len(hinge_losses) > 0:\n    ax.hist(hinge_losses, bins=30, color='coral', edgecolor='white', alpha=0.8)\nax.axvline(avg_hinge, color='red', linestyle='--', linewidth=2, label=f'Mean={avg_hinge:.3f}')\nax.set_xlabel('Hinge Loss'); ax.set_title('Hinge Loss Distribution', fontweight='bold'); ax.legend()\n\n# (f) Root vs Part Contribution\nax = fig.add_subplot(gs[1, 2])\nax.pie([root_contrib_pct, part_contrib_pct],\n       labels=[f'Root Filter\\n{root_contrib_pct:.1f}%', f'Part Filters\\n{part_contrib_pct:.1f}%'],\n       colors=['#e74c3c', '#3498db'], autopct='%1.1f%%', startangle=90,\n       wedgeprops={'edgecolor': 'white', 'linewidth': 2})\nax.set_title('Root vs Part Contribution\\n(Weight Norm Ratio)', fontweight='bold')\n\n# (g) Hypothesis Score Distribution\nax = fig.add_subplot(gs[2, 0])\nif len(scores_arr) > 0:\n    if np.any(y_val_binary == 1):\n        ax.hist(root_scores_pos, bins=20, alpha=0.6, color='green', label='TP (IoU≥0.5)', edgecolor='white')\n    if np.any(y_val_binary == -1):\n        ax.hist(root_scores_neg, bins=20, alpha=0.6, color='red', label='FP (IoU<0.5)', edgecolor='white')\nax.set_xlabel('SVM Score'); ax.set_title('Hypothesis Score Distribution', fontweight='bold'); ax.legend()\n\n# (h) Deformation Cost Distribution\nax = fig.add_subplot(gs[2, 1])\nax.hist(deform_costs, bins=20, color='purple', edgecolor='white', alpha=0.8)\nax.axvline(avg_deform, color='red', linestyle='--', label=f'Mean={avg_deform:.4f}')\nax.set_xlabel('Deformation Penalty'); ax.set_title('Deformation Cost Distribution', fontweight='bold'); ax.legend()\n\n# (i) Summary table\nax = fig.add_subplot(gs[2, 2])\nax.axis('off')\ntable_data = [\n    ['Average Precision (AP@0.5)', f'{ap_05:.4f}'],\n    ['Average Precision (AP@0.3)', f'{ap_03:.4f}'],\n    ['Precision', f'{precision_val:.4f}'],\n    ['Recall', f'{recall_val:.4f}'],\n    ['F1 Score', f'{f1_val:.4f}'],\n    ['F2 Score', f'{f2_val:.4f}'],\n    ['Mean IoU', f'{avg_iou:.4f}'],\n    ['Avg Hinge Loss', f'{avg_hinge:.4f}'],\n    ['Root Filter Score (mean TP)', f'{np.mean(root_scores_pos):.4f}'],\n    ['Root Contribution', f'{root_contrib_pct:.1f}%'],\n    ['Part Contribution', f'{part_contrib_pct:.1f}%'],\n    ['Avg Deformation Cost', f'{avg_deform:.4f}'],\n]\ntable = ax.table(cellText=table_data, colLabels=['Metric', 'Value'], loc='center', cellLoc='left')\ntable.auto_set_font_size(False); table.set_fontsize(9); table.scale(1, 1.4)\nfor (row, col), cell in table.get_celld().items():\n    if row == 0:\n        cell.set_facecolor('#2c3e50'); cell.set_text_props(color='white', fontweight='bold')\n    else:\n        cell.set_facecolor('#ecf0f1' if row % 2 == 0 else 'white')\nax.set_title('Metrics Summary', fontweight='bold', pad=20)\n\nplt.suptitle(f'Phase 7 — Metrics Dashboard ({TARGET_CLASS})', fontsize=16, fontweight='bold', y=1.01)\nplt.savefig(OUTPUT_DIR / 'metrics_dashboard.png', dpi=150, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"2090b26e-8e51-4d19-a83f-766320d5d120","cell_type":"markdown","source":"### 📊 Visualization 10 — DPM Interpretability (Root Response Heatmap + Parts)","metadata":{}},{"id":"8e06f83d-aca4-4456-af59-a64cc2bb21fc","cell_type":"code","source":"def visualise_dpm_detection(image_id, detections, part_anchors, class_name,\n                             root_svm, feat_scaler, save=True):\n    \"\"\"4-panel interpretability visualisation with root heatmap + parts.\"\"\"\n    img = load_dicom(image_id)\n    H, W = img.shape\n    _, hog_vis = HOGExtractor().extract(img)\n\n    # Sparse response map\n    ps, stride = 64, 32\n    nh = (H - ps) // stride + 1\n    nw = (W - ps) // stride + 1\n    resp = np.zeros((nh, nw), dtype=np.float32)\n    for ri in range(nh):\n        for ci in range(nw):\n            y0, x0 = ri*stride, ci*stride\n            patch   = img[y0:y0+ps, x0:x0+ps]\n            f, _    = hog(patch, orientations=9, pixels_per_cell=(8,8),\n                          cells_per_block=(2,2), visualize=True, channel_axis=None)\n            resp[ri, ci] = feat_scaler.transform([f])[0] @ root_svm.coef_[0]\n\n    fig, axes = plt.subplots(1, 4, figsize=(22, 5))\n    axes[0].imshow(img, cmap='gray'); axes[0].set_title('Original CXR'); axes[0].axis('off')\n    axes[1].imshow(hog_vis, cmap='hot'); axes[1].set_title('HOG Features'); axes[1].axis('off')\n\n    axes[2].imshow(img, cmap='gray', alpha=0.6)\n    axes[2].imshow(cv2.resize(resp, (W, H)), cmap='RdYlGn', alpha=0.5)\n    axes[2].set_title('Root Filter Response\\n(Heatmap)'); axes[2].axis('off')\n\n    axes[3].imshow(img, cmap='gray')\n    colors = ['red', 'cyan', 'lime', 'orange', 'magenta']\n    for i, (score, x1, y1, x2, y2) in enumerate(detections[:1]):\n        rect = plt.Rectangle((x1, y1), x2-x1, y2-y1,\n                              linewidth=3, edgecolor='red', facecolor='none')\n        axes[3].add_patch(rect)\n        axes[3].text(x1, y1-8, f'ROOT: {score:.2f}', color='red', fontsize=8,\n                     bbox=dict(boxstyle='round', facecolor='white', alpha=0.7))\n        for p_idx, (rx, ry, fw, fh) in enumerate(part_anchors):\n            px1 = int(x1 + rx*(x2-x1)); py1 = int(y1 + ry*(y2-y1))\n            px2 = int(px1 + fw*(x2-x1)); py2 = int(py1 + fh*(y2-y1))\n            prect = plt.Rectangle((px1, py1), px2-px1, py2-py1,\n                                   linewidth=2, edgecolor=colors[p_idx % len(colors)],\n                                   facecolor=colors[p_idx % len(colors)],\n                                   linestyle='--', alpha=0.12)\n            axes[3].add_patch(prect)\n            prect2 = plt.Rectangle((px1, py1), px2-px1, py2-py1,\n                                    linewidth=2, edgecolor=colors[p_idx % len(colors)],\n                                    facecolor='none', linestyle='--')\n            axes[3].add_patch(prect2)\n            axes[3].text(px1, py1-3, f'P{p_idx}', color=colors[p_idx % len(colors)], fontsize=7,\n                         fontweight='bold')\n    axes[3].set_title('DPM: Root + Parts'); axes[3].axis('off')\n\n    plt.suptitle(f'Interpretable DPM Detection — {class_name}', fontsize=13, fontweight='bold')\n    plt.tight_layout()\n    if save:\n        plt.savefig(OUTPUT_DIR / f'dpm_{image_id}_{class_name.replace(\" \",\"_\")}.png', dpi=150)\n    plt.show()\n\n\n# ── Run visualisation ──\nif dets:\n    visualise_dpm_detection(\n        val_row['image_id'], dets, PART_ANCHORS[TARGET_CLASS],\n        TARGET_CLASS, saved_model['svm'], saved_model['scaler']\n    )\nelse:\n    print('No detections found above threshold for this image.')","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"ad0aad43-670c-431d-9a67-04e5d15db19a","cell_type":"markdown","source":"## 🔍 Anomaly Detection Results — Final Summary\n","metadata":{}},{"id":"cf4d9935-62bc-49a4-a1e2-b1912afc1cfc","cell_type":"code","source":"print('='*70)\nprint('         ANOMALY DETECTION RESULTS SUMMARY')\nprint('='*70)\nprint(f'\\nTarget Anomaly Class: {TARGET_CLASS}')\nprint(f'Images Evaluated: {n_images_evaluated}')\nprint(f'Total Ground Truth Annotations: {n_gt_total}')\nprint(f'Total Detections: {len(all_preds)}')\nprint(f'True Positives: {tp_total}')\nprint(f'False Positives: {fp_total}')\nprint(f'False Negatives: {fn_total}')\nprint()\nprint('─' * 50)\nprint('  ANOMALIES FOUND:')\nprint('─' * 50)\n\nif tp_total > 0:\n    matched_imgs = set()\n    sorted_p = sorted(all_preds, key=lambda x: -x[0])\n    for score, x1, y1, x2, y2, img_id in sorted_p:\n        best_iou = max([compute_iou((x1,y1,x2,y2), gb) for gb in gt_dict.get(img_id, [])] or [0])\n        if best_iou >= 0.5 and img_id not in matched_imgs:\n            print(f'  ✅ {TARGET_CLASS} DETECTED in image {img_id}')\n            print(f'     Location: ({x1},{y1}) → ({x2},{y2})  |  Score: {score:.3f}  |  IoU: {best_iou:.3f}')\n            matched_imgs.add(img_id)\n    print(f'\\n  Total confirmed anomalies: {len(matched_imgs)} images')\nelse:\n    print('  ⚠️  No confirmed anomalies detected at IoU≥0.5 threshold')\n\nprint()\nprint('─' * 50)\nprint('  PERFORMANCE METRICS:')\nprint('─' * 50)\nprint(f'  Average Precision (AP@0.5):  {ap_05:.4f}')\nprint(f'  Average Precision (AP@0.3):  {ap_03:.4f}')\nprint(f'  Precision:                   {precision_val:.4f}')\nprint(f'  Recall:                      {recall_val:.4f}')\nprint(f'  F1 Score:                    {f1_val:.4f}')\nprint(f'  F2 Score:                    {f2_val:.4f}')\nprint(f'  Mean BBox IoU:               {avg_iou:.4f}')\nprint(f'  Avg Hinge Loss:              {avg_hinge:.4f}')\nprint(f'  Root Filter Score (TP mean):  {np.mean(root_scores_pos):.4f}')\nprint(f'  Part Contribution:           {part_contrib_pct:.1f}%')\nprint(f'  Avg Deformation Cost:        {avg_deform:.4f}')\nprint()\nprint('  FROC Sensitivities:')\nfor fpi, sens in zip(fp_per_img_rates, froc_sensitivities):\n    print(f'    @ {fpi:.1f} FP/image: {sens:.3f}')\nprint('='*70)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"66bcd9fd-8c7d-4eb3-af6b-54a328d5a9de","cell_type":"markdown","source":"## Phase 8 — Kaggle Submission\n","metadata":{}},{"id":"82dd8bd7-c721-4f8c-b696-4dd076993b45","cell_type":"code","source":"test_csv    = pd.read_csv(TEST_CSV)\ntest_img_ids = test_csv['image_id'].unique()\nprint(f'Test images: {len(test_img_ids)}')","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"55962a05-754b-4859-bcc1-3a2d11bee71d","cell_type":"code","source":"submission_rows = []\n\nfor image_id in tqdm(test_img_ids, desc='Generating submission'):\n    found_any = False\n    for cls in TOP_CLASSES:\n        model_path = MODEL_DIR / f\"root_filter_{cls.replace(' ', '_')}.pkl\"\n        if not model_path.exists():\n            continue\n        m    = joblib.load(model_path)\n        img  = load_dicom(image_id, split='test')\n        dets = sliding_window_detect(img, m['svm'], m['scaler'], stride=32, threshold=0.3)\n        dets = nms(dets)\n        for score, x1, y1, x2, y2 in dets[:3]:\n            submission_rows.append({\n                'image_id':   image_id,\n                'class_name': cls,\n                'score':      round(float(score), 4),\n                'x_min': x1, 'y_min': y1, 'x_max': x2, 'y_max': y2\n            })\n            found_any = True\n    if not found_any:\n        submission_rows.append({\n            'image_id':   image_id,\n            'class_name': 'No finding',\n            'score':      1.0,\n            'x_min': 0, 'y_min': 0, 'x_max': 1, 'y_max': 1\n        })\n\nsubmission_df = pd.DataFrame(submission_rows)\nsubmission_df.to_csv('/kaggle/working/submission.csv', index=False)\nprint('Submission saved!')\nprint(submission_df['class_name'].value_counts())\nprint(submission_df.head())","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"f5b60aa3-2056-47bb-958e-036b10a8ad7d","cell_type":"markdown","source":"## ✅ Summary\n\n| Phase | Status |\n|---|---|\n| 0 — Environment Setup | ✅ |\n| 1 — Data Preparation & EDA | ✅ + Visualization |\n| 2 — HOG + Image Pyramid | ✅ + Visualization |\n| 3 — Root Filter (LinearSVC) | ✅ + Visualization |\n| 4 — Part Filters + Latent SVM | ✅ + Root/Part Decomposition |\n| 5 — Deformation Model | ✅ + Cost Surface Plot |\n| 6 — Sliding Window + NMS | ✅ + Before/After NMS |\n| 7 — Evaluation + Metrics Dashboard | ✅ All 11 metrics |\n| 8 — Kaggle Submission | ✅ |\n\n### Metrics Computed:\n- Average Precision (AP@0.5, AP@0.3)\n- Precision & Recall\n- F1 Score & F2 Score\n- Bounding Box IoU Distribution\n- Hypothesis Score (SVM decision function)\n- Deformation Cost Distribution\n- Hinge Loss\n- FROC Curve\n- Root Filter Score\n- Part Contribution Score\n","metadata":{}}]}