{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","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":"gpu","dataSources":[{"sourceType":"competition","sourceId":10338,"databundleVersionId":862042}],"dockerImageVersionId":31287,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"!pip install faiss-cpu","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-24T07:43:10.614928Z","iopub.execute_input":"2026-02-24T07:43:10.615206Z","iopub.status.idle":"2026-02-24T07:43:15.871986Z","shell.execute_reply.started":"2026-02-24T07:43:10.615182Z","shell.execute_reply":"2026-02-24T07:43:15.871310Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nimport cv2\nimport torch\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom torch.utils.data import Dataset, DataLoader\nimport torchvision.transforms as transforms\nfrom sklearn.neighbors import NearestNeighbors\nimport torch.nn.functional as F","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-24T07:43:15.873682Z","iopub.execute_input":"2026-02-24T07:43:15.873918Z","iopub.status.idle":"2026-02-24T07:43:24.580785Z","shell.execute_reply.started":"2026-02-24T07:43:15.873893Z","shell.execute_reply":"2026-02-24T07:43:24.580172Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Configuration\nDEVICE = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\nBATCH_SIZE = 2\nIMG_SIZE = 112\nTOP_K = 10 # For Top-K pooling\nNUM_WORKERS = 1\nPATCH_SIZE = 14","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-24T07:43:24.581559Z","iopub.execute_input":"2026-02-24T07:43:24.581909Z","iopub.status.idle":"2026-02-24T07:43:24.640959Z","shell.execute_reply.started":"2026-02-24T07:43:24.581887Z","shell.execute_reply":"2026-02-24T07:43:24.640362Z"},"_kg_hide-input":false},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import torch, numpy as np, random\n\ntorch.manual_seed(42)\nnp.random.seed(42)\nrandom.seed(42)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-24T07:43:24.642536Z","iopub.execute_input":"2026-02-24T07:43:24.642896Z","iopub.status.idle":"2026-02-24T07:43:24.664609Z","shell.execute_reply.started":"2026-02-24T07:43:24.642868Z","shell.execute_reply":"2026-02-24T07:43:24.664086Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nos.environ[\"PYTORCH_CUDA_ALLOC_CONF\"] = \"max_split_size_mb:64\"\n\n# The presentation specifies converting grayscale to 3-channel and normalizing.\n\ntransform = transforms.Compose([\n    transforms.ToPILImage(),\n    transforms.Resize((IMG_SIZE, IMG_SIZE)),\n    transforms.Grayscale(num_output_channels=3), # Intensity replication\n    transforms.ToTensor(),\n    transforms.Normalize(mean=[0.485, 0.456, 0.406], std=[0.229, 0.224, 0.225]) #\n])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-24T07:43:24.665509Z","iopub.execute_input":"2026-02-24T07:43:24.665740Z","iopub.status.idle":"2026-02-24T07:43:24.670046Z","shell.execute_reply.started":"2026-02-24T07:43:24.665719Z","shell.execute_reply":"2026-02-24T07:43:24.669460Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class RSNADataset(Dataset):\n    def __init__(self, df, img_dir, transform=None, is_train=True):\n        self.img_dir = img_dir\n        self.transform = transform\n        \n        if is_train:\n            self.df = df[df['class'] == 'Normal'].reset_index(drop=True)\n        else:\n            self.df = df.reset_index(drop=True)\n            \n    def __len__(self):\n        return len(self.df)\n\n    def __getitem__(self, idx):\n        img_name = os.path.join(self.img_dir, f\"{self.df.loc[idx, 'patientId']}.dcm\")\n        \n        import pydicom\n        if not os.path.exists(img_name):\n            raise FileNotFoundError(img_name)\n\n        image = pydicom.dcmread(img_name).pixel_array\n        if image.ndim > 2:\n            image = image[..., 0]\n\n        image = cv2.normalize(image, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)\n\n        if self.transform:\n            image = self.transform(image)\n                \n        label = 0 if self.df.loc[idx, 'class'] == 'Normal' else 1\n        return image, label","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-24T07:43:24.670938Z","iopub.execute_input":"2026-02-24T07:43:24.671214Z","iopub.status.idle":"2026-02-24T07:43:24.682274Z","shell.execute_reply.started":"2026-02-24T07:43:24.671182Z","shell.execute_reply":"2026-02-24T07:43:24.681624Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Initializing DINOv2 (ViT-B/14) backbone.\n\nimport timm\n\nprint(\"Loading DINOv2 ViT-B/14 model (timm)\")\ndinov2_model = timm.create_model(\n    \"vit_base_patch14_dinov2.lvd142m\",\n    pretrained = True,\n    img_size = IMG_SIZE\n)\n\ndinov2_model = dinov2_model.to(DEVICE)\ndinov2_model.eval()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-24T07:43:24.683084Z","iopub.execute_input":"2026-02-24T07:43:24.683360Z","iopub.status.idle":"2026-02-24T07:43:35.730865Z","shell.execute_reply.started":"2026-02-24T07:43:24.683339Z","shell.execute_reply":"2026-02-24T07:43:35.730141Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Extract patch-level features to capture long-range dependencies.\n\nfrom tqdm import tqdm\nimport os\n\ndef extract_patch_features_to_disk(dataloader, model, save_path, desc=\"Extracting\"):\n    model.eval()\n    all_files = []\n    \n    with torch.no_grad():\n        for i, (images, _) in enumerate(tqdm(dataloader, desc=desc)):\n            images = images.to(DEVICE)\n            feats = model.forward_features(images)\n\n            if isinstance(feats, dict):\n                patch_tokens = feats[\"x_norm_patchtokens\"]   # timm dict output\n            else:\n                patch_tokens = feats[:, 1:, :]               # tensor output (remove CLS token)\n            \n            # REMOVED the mean(dim=1) pooling here. \n            # Keep the shape as (Batch, Num_Patches, 768)\n            batch_feats = patch_tokens.cpu().numpy()\n            \n            fname = f\"{save_path}_batch_{i}.npy\"\n            np.save(fname, batch_feats)\n            all_files.append(fname)\n\n            # free memory aggressively\n            del images, feats, patch_tokens, batch_feats\n            torch.cuda.empty_cache()\n\n    return all_files","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-24T07:43:35.731989Z","iopub.execute_input":"2026-02-24T07:43:35.732460Z","iopub.status.idle":"2026-02-24T07:43:35.738134Z","shell.execute_reply.started":"2026-02-24T07:43:35.732437Z","shell.execute_reply":"2026-02-24T07:43:35.737568Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Measure distance between test features and the learned healthy distribution.\n\nimport faiss\nfrom tqdm import tqdm\nimport numpy as np\n\ndef compute_anomaly_scores_faiss(train_features, test_features, batch_size=128):\n    \"\"\"\n    Fast FAISS CPU k-NN using IVF (approximate search).\n    Cosine similarity via normalized Inner Product.\n    \"\"\"\n\n    print(\"Preparing FAISS IVF index...\")\n\n    # Flatten train features\n    num_train_imgs, num_train_patches, feat_dim = train_features.shape\n    train_flat = train_features.reshape(-1, feat_dim).astype(np.float32)\n\n    # Normalize for cosine similarity\n    faiss.normalize_L2(train_flat)\n\n    # -----------------------\n    # IVF PARAMETERS\n    # -----------------------\n    nlist = 4096  # number of clusters (try 2048 if memory tight)\n\n    quantizer = faiss.IndexFlatIP(feat_dim)\n    index = faiss.IndexIVFFlat(\n        quantizer,\n        feat_dim,\n        nlist,\n        faiss.METRIC_INNER_PRODUCT\n    )\n\n    print(\"Training IVF index (one-time cost)...\")\n    index.train(train_flat)\n\n    print(\"Adding vectors to index...\")\n    index.add(train_flat)\n\n    # Controls speed/accuracy tradeoff\n    index.nprobe = 16  # try 8–32\n\n    print(\"Index size:\", index.ntotal)\n    print(\"Running FAISS IVF k-NN search...\")\n\n    num_test_imgs, num_test_patches, _ = test_features.shape\n    all_scores = []\n\n    for i in tqdm(range(0, num_test_imgs, batch_size), desc=\"FAISS IVF k-NN\"):\n        batch_test = test_features[i:i + batch_size]\n        B, P, D = batch_test.shape\n\n        batch_flat = batch_test.reshape(-1, D).astype(np.float32)\n        faiss.normalize_L2(batch_flat)\n\n        sims, _ = index.search(batch_flat, 1)\n\n        dists = 1.0 - sims\n        batch_scores = dists.reshape(B, P)\n        all_scores.append(batch_scores)\n\n    return np.concatenate(all_scores, axis=0)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-24T07:43:35.739014Z","iopub.execute_input":"2026-02-24T07:43:35.739294Z","iopub.status.idle":"2026-02-24T07:43:35.787126Z","shell.execute_reply.started":"2026-02-24T07:43:35.739263Z","shell.execute_reply":"2026-02-24T07:43:35.786612Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def load_features_from_disk(file_list):\n    feats = [np.load(f) for f in file_list]\n    return np.concatenate(feats, axis=0)\n\ndef classify_and_localize(patch_scores, img_size=IMG_SIZE, patch_size=14):\n    num_imgs, num_patches = patch_scores.shape\n    grid_size = img_size // patch_size\n    \n    # --- Classification ---\n    # Aggregate patch-level scores using top-k pooling\n    sorted_scores = np.sort(patch_scores, axis=1)[:, ::-1] # Sort descending\n    global_anomaly_scores = np.mean(sorted_scores[:, :TOP_K], axis=1) # Top-K average\n    \n    # --- Localization ---\n    # Project patch-level scores back to original image resolution\n    heatmaps = []\n    for score in patch_scores:\n        # Reshape to 2D grid\n        spatial_score = score.reshape(grid_size, grid_size)\n        \n        # Resize to original image size for pixel-level heatmap\n        heatmap = cv2.resize(spatial_score, (img_size, img_size), interpolation=cv2.INTER_CUBIC)\n        \n        # Normalize heatmap for visualization\n        heatmap = (heatmap - heatmap.min()) / (heatmap.max() - heatmap.min() + 1e-8)\n        heatmaps.append(heatmap)\n        \n    return global_anomaly_scores, np.array(heatmaps)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-24T07:43:35.788712Z","iopub.execute_input":"2026-02-24T07:43:35.788930Z","iopub.status.idle":"2026-02-24T07:43:35.794455Z","shell.execute_reply.started":"2026-02-24T07:43:35.788910Z","shell.execute_reply":"2026-02-24T07:43:35.793856Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from sklearn.model_selection import train_test_split\nfrom torch.utils.data import DataLoader\n\nif __name__ == \"__main__\":\n\n    # 1. Load dataframe\n    df = pd.read_csv(\n        '/kaggle/input/rsna-pneumonia-detection-challenge/stage_2_train_labels.csv'\n    )\n\n    df['class'] = df['Target'].map({0: 'Normal', 1: 'Pneumonia'})\n    img_dir = '/kaggle/input/rsna-pneumonia-detection-challenge/stage_2_train_images'\n\n    train_df, test_df = train_test_split(\n        df,\n        test_size=0.2,\n        stratify=df['class'],\n        random_state=42\n    )\n\n    print(\"Train size:\", len(train_df))\n    print(\"Test size:\", len(test_df))\n\n    # Train dataset → ONLY Normal images from TRAIN split\n    train_dataset = RSNADataset(\n        train_df,\n        img_dir,\n        transform=transform,\n        is_train=True\n    )\n\n    # Test dataset → ALL images from TEST split\n    test_dataset = RSNADataset(\n        test_df,\n        img_dir,\n        transform=transform,\n        is_train=False\n    )\n\n    train_loader = DataLoader(\n        train_dataset,\n        batch_size=BATCH_SIZE,\n        shuffle=True,\n        num_workers=NUM_WORKERS\n    )\n\n    test_loader = DataLoader(\n        test_dataset,\n        batch_size=BATCH_SIZE,\n        shuffle=False,\n        num_workers=NUM_WORKERS\n    )\n    \n    import time\n\n    print(\"Extracting healthy (train) features...\")\n    start = time.time()\n\n    train_files = extract_patch_features_to_disk(\n        train_loader,\n        dinov2_model,\n        \"train_feats\",\n        \"Train features\"\n    )\n\n    train_features = load_features_from_disk(train_files)\n\n    print(\"Train feature time (min):\",\n          (time.time() - start) / 60)\n\n    print(\"Extracting test features...\")\n    test_files = extract_patch_features_to_disk(\n        test_loader,\n        dinov2_model,\n        \"test_feats\",\n        \"Test features\"\n    )\n\n    test_features = load_features_from_disk(test_files)\n\n    patch_scores = compute_anomaly_scores_faiss(\n        train_features,\n        test_features,\n        batch_size=64\n    )\n\n    print(\"Patch scores shape:\", patch_scores.shape)\n\n    image_scores, heatmaps = classify_and_localize(\n        patch_scores,\n        img_size=IMG_SIZE,\n        patch_size=PATCH_SIZE\n    )\n\n    print(\"Pipeline finished successfully.\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-24T07:43:35.795285Z","iopub.execute_input":"2026-02-24T07:43:35.795505Z","iopub.status.idle":"2026-02-24T08:03:37.469212Z","shell.execute_reply.started":"2026-02-24T07:43:35.795481Z","shell.execute_reply":"2026-02-24T08:03:37.468359Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nfrom sklearn.metrics import roc_auc_score, roc_curve, accuracy_score, f1_score, confusion_matrix\nimport matplotlib.pyplot as plt\n\n# 1. Instantly grab labels directly from the dataframe (no disk I/O)\ny_true = (test_dataset.df['class'] != 'Normal').astype(int).values\n\n# 2. Calculate AUC \nauc = roc_auc_score(y_true, image_scores)\nprint(f\"AUC (%): {auc * 100:.2f}\")\n\n# \n# 3. Youden’s J statistic (optimal ROC threshold)\nfpr, tpr, thresholds = roc_curve(y_true, image_scores)\nj_scores = tpr - fpr\nbest_thresh = thresholds[np.argmax(j_scores)]\n\ny_pred = (image_scores > best_thresh).astype(int)\nprint(f\"Best threshold: {best_thresh:.4f}\")\n\n# 4. Standard Metrics\nacc = accuracy_score(y_true, y_pred)\nf1  = f1_score(y_true, y_pred)\n\nprint(f\"Accuracy (%): {acc * 100:.2f}\")\nprint(f\"F1-score (%): {f1 * 100:.2f}\")\n\ncm = confusion_matrix(y_true, y_pred)\nprint(\"Confusion Matrix:\\n\", cm)\n\ndef visualize_anomalies(dataset, heatmaps, y_true, y_pred, num_samples=3):\n    \"\"\"Plots original images alongside their anomaly heatmaps[cite: 74, 133].\"\"\"\n    # Find a few True Positives (Actual Anomaly, Predicted Anomaly)\n    tp_indices = np.where((y_true == 1) & (y_pred == 1))[0]\n    \n    if len(tp_indices) == 0:\n        print(\"No True Positives found to visualize.\")\n        return\n        \n    samples = np.random.choice(tp_indices, min(num_samples, len(tp_indices)), replace=False)\n    \n    fig, axes = plt.subplots(len(samples), 2, figsize=(10, 5 * len(samples)))\n    if len(samples) == 1: axes = [axes]\n    \n    for idx, ax_row in zip(samples, axes):\n        # We can safely use dataset[idx] here because we only load a few images\n        img_tensor, label = dataset[idx] \n        \n        # Convert tensor back to numpy for plotting\n        # Un-normalize the image for viewing\n        img_np = img_tensor.permute(1, 2, 0).numpy()\n        mean = np.array([0.485, 0.456, 0.406])\n        std = np.array([0.229, 0.224, 0.225])\n        img_np = std * img_np + mean\n        img_np = np.clip(img_np, 0, 1)\n        \n        heatmap = heatmaps[idx]\n        \n        ax_row[0].imshow(img_np)\n        ax_row[0].set_title(f\"Original X-Ray (Index {idx})\")\n        ax_row[0].axis('off')\n        \n        # Overlay heatmap\n        ax_row[1].imshow(img_np)\n        im = ax_row[1].imshow(heatmap, cmap='jet', alpha=0.5) # Alpha blends the heatmap \n        ax_row[1].set_title(\"Anomaly Heatmap Overlay\")\n        ax_row[1].axis('off')\n        fig.colorbar(im, ax=ax_row[1], fraction=0.046, pad=0.04)\n        \n    plt.tight_layout()\n    plt.show()\n\n# Run the visualization\nprint(\"\\nGenerating heatmaps for True Positives...\")\nvisualize_anomalies(test_dataset, heatmaps, y_true, y_pred)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-24T08:03:37.470814Z","iopub.execute_input":"2026-02-24T08:03:37.471412Z","iopub.status.idle":"2026-02-24T08:03:38.225954Z","shell.execute_reply.started":"2026-02-24T08:03:37.471382Z","shell.execute_reply":"2026-02-24T08:03:38.225149Z"}},"outputs":[],"execution_count":null}]}