{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":20561,"databundleVersionId":1187889,"sourceType":"competition"},{"sourceId":13060062,"sourceType":"datasetVersion","datasetId":8270323},{"sourceId":13104418,"sourceType":"datasetVersion","datasetId":8300843},{"sourceId":262712410,"sourceType":"kernelVersion"}],"dockerImageVersionId":31089,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Covid-19 Segmentation\n\nJesus Ramirez Delgado A01274723\n\nGrant Nathaniel Keegan A01700753\n\nLuis Adrian Uribe Cruz A01783129\n\nFidel Alexander Bonilla Montalvo A01798199","metadata":{}},{"cell_type":"markdown","source":"## 1. Introduction\n\nThe COVID-19 pandemic posed one of the most significant recent challenges to healthcare systems worldwide. One of the most widely used methods for diagnosing and monitoring patients was chest computed tomography (CT), which enables the identification of disease-specific patterns. However, manually interpreting hundreds of tomographic slices is a slow, radiologist-dependent, and variable process, opening the door to automatic decision-support tools.\n\nIn this context, the objective was to develop machine-learning models capable of automatically segmenting COVID-19–affected regions in axial CT slices. The core task is semantic segmentation—that is, classifying every pixel in the image as belonging to a specific clinical class. The challenge is to evaluate the models’ ability to accurately identify these lesions while optimizing metrics such as the Dice coefficient and Intersection over Union (IoU), which are especially relevant under severe class imbalance.\n\nThis task is relevant in two ways: on the one hand, it contributes a potentially useful tool for clinical practice, helping physicians quickly detect critical patterns in patients with COVID-19; on the other, it strengthens research in medical computer vision by providing standardized datasets and fostering open comparisons of methodologies.\n\nTo address this problem, we used a reproducible, Jupyter-notebook-based workflow in which we developed an end-to-end pipeline that includes:\n\nExploratory and preliminary analysis of the datasets.\n\nAn ETL process that normalizes images, selects relevant classes, and prevents information leakage when splitting the data.\n\nImplementation and training of a 2D U-Net model in PyTorch, tuned for medical segmentation with limited data.\n\nVisualization and validation of results, both quantitative and qualitative.","metadata":{}},{"cell_type":"markdown","source":"## 2. Data Analysis - Grant ","metadata":{}},{"cell_type":"markdown","source":"### 2.1 Exploratory Data Analysis (EDA)","metadata":{}},{"cell_type":"markdown","source":"The dataset is made up of 4 NumPy arrays (,npy) to train the model:\n- images_medseg.npy\n- images_radiopedia.npy\n- masks_medseg.npy\n- masks_radiopedia.npy\n\nAnd one for testing:\n- test_images_medseg.npy\n\nThe data sets can be separated by sorce:\n- medseg\n- radiopedia\n\nCategories\n\nEach source is split into two categories:\n\n- images: axial CT slices stored as grayscale intensities in Hounsfield Units (HU).\nTypical shape: (N, H, W, 1) or (N, H, W); dtype usually float32/float64.\nCommon lung window for visualization: WL≈-600, WW≈1500.\n\n- masks: per-pixel segmentation labels aligned slice-by-slice with images from the same source.\nFormat can be either:\n\n    - Label map: (N, H, W) with integer classes {0,…,K-1}, or\n\n    - One-hot: (N, H, W, K) where channel K indexes classes.\n      In MedSeg, the canonical classes are:\n\n    0: ground glass\n    \n    1: consolidation\n    \n    2: lungs other\n    \n    3: background","metadata":{}},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session\n\n# Bandera global\nRUN_TRAIN = False  # pon True cuando sí quieras entrenar\n\n# Magic para saltar celdas condicionalmente\nfrom IPython.core.magic import register_cell_magic\n\n@register_cell_magic\ndef skip_if(line, cell):\n    if eval(line, globals(), locals()):\n        print(f\"Skiped block: {line}\")\n    else:\n        exec(cell, globals(), locals())\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"prefix = '/kaggle/input/covid-segmentation/'\n\nimages_radiopedia = np.load(os.path.join(prefix, 'images_radiopedia.npy')).astype(np.float32)\nmasks_radiopedia = np.load(os.path.join(prefix, 'masks_radiopedia.npy')).astype(np.int8)\nimages_medseg = np.load(os.path.join(prefix, 'images_medseg.npy')).astype(np.float32)\nmasks_medseg = np.load(os.path.join(prefix, 'masks_medseg.npy')).astype(np.int8)\n\ntest_images_medseg = np.load(os.path.join(prefix, 'test_images_medseg.npy')).astype(np.float32)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-18T12:43:48.226158Z","iopub.execute_input":"2025-09-18T12:43:48.226386Z","iopub.status.idle":"2025-09-18T12:43:50.039782Z","shell.execute_reply.started":"2025-09-18T12:43:48.226368Z","shell.execute_reply":"2025-09-18T12:43:50.038999Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\n\ndef shape_NHWK(arr, is_image: bool):\n    \"\"\"Return (N,H,W,K) for images/masks from common shapes.\"\"\"\n    if arr.ndim == 4:\n        N, H, W, K = arr.shape\n        return N, H, W, K\n    if arr.ndim == 3:\n        N, H, W = arr.shape\n        K = 1 if is_image else 1\n        return N, H, W, K\n    raise ValueError(f\"Unexpected shape {arr.shape} (is_image={is_image})\")\n\ndef describe_K(array_kind: str, K: int):\n    if array_kind == \"Images\":\n        return f\"{K} (grayscale image)\" if K == 1 else str(K)\n    # Masks\n    if K == 4:\n        return \"4 (0=GGO, 1=Consolidation, 2=Lungs-other, 3=Background)\"\n    if K == 2:\n        return \"2 (0=GGO, 1=Consolidation)\"\n    if K == 1:\n        return \"1 (single-label mask)\"\n    return str(K)\n\nrows = []\n\n# MedSeg\nN,H,W,K = shape_NHWK(images_medseg, is_image=True)\nrows.append({\n    \"Split / Source\": \"MedSeg\",\n    \"Array\": \"Images\",\n    \"Shape (N, H, W, K)\": f\"({N}, {H}, {W}, {K})\",\n    \"N\": f\"{N} (slices)\",\n    \"H×W\": f\"{H}×{W}\",\n    \"K (channels/classes)\": describe_K(\"Images\", K),\n})\n\nN,H,W,K = shape_NHWK(masks_medseg, is_image=False)\nrows.append({\n    \"Split / Source\": \"MedSeg\",\n    \"Array\": \"Masks\",\n    \"Shape (N, H, W, K)\": f\"({N}, {H}, {W}, {K})\",\n    \"N\": f\"{N} (slices)\",\n    \"H×W\": f\"{H}×{W}\",\n    \"K (channels/classes)\": describe_K(\"Masks\", K),\n})\n\n# MedSeg (test)\nN,H,W,K = shape_NHWK(test_images_medseg, is_image=True)\nrows.append({\n    \"Split / Source\": \"MedSeg (Test)\",\n    \"Array\": \"Images\",\n    \"Shape (N, H, W, K)\": f\"({N}, {H}, {W}, {K})\",\n    \"N\": f\"{N} (slices)\",\n    \"H×W\": f\"{H}×{W}\",\n    \"K (channels/classes)\": describe_K(\"Images\", K),\n})\n\n# Radiopaedia\nN,H,W,K = shape_NHWK(images_radiopedia, is_image=True)\nrows.append({\n    \"Split / Source\": \"Radiopaedia\",\n    \"Array\": \"Images\",\n    \"Shape (N, H, W, K)\": f\"({N}, {H}, {W}, {K})\",\n    \"N\": f\"{N} (slices)\",\n    \"H×W\": f\"{H}×{W}\",\n    \"K (channels/classes)\": describe_K(\"Images\", K),\n})\n\nN,H,W,K = shape_NHWK(masks_radiopedia, is_image=False)\nrows.append({\n    \"Split / Source\": \"Radiopaedia\",\n    \"Array\": \"Masks\",\n    \"Shape (N, H, W, K)\": f\"({N}, {H}, {W}, {K})\",\n    \"N\": f\"{N} (slices)\",\n    \"H×W\": f\"{H}×{W}\",\n    \"K (channels/classes)\": describe_K(\"Masks\", K),\n})\n\n# Optional note row\nrows.append({\n    \"Split / Source\": \"Notes\",\n    \"Array\": \"Meaning of dimensions\",\n    \"Shape (N, H, W, K)\": \"\",\n    \"N\": \"Number of slices\",\n    \"H×W\": \"Height × Width (pixels)\",\n    \"K (channels/classes)\": \"Images: 1 channel; Masks: 4 classes (as listed)\",\n})\n\ndf = pd.DataFrame(rows, columns=[\n    \"Split / Source\", \"Array\", \"Shape (N, H, W, K)\", \"N\", \"H×W\", \"K (channels/classes)\"\n])\n\ndisplay(df)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-18T12:43:50.040627Z","iopub.execute_input":"2025-09-18T12:43:50.040885Z","iopub.status.idle":"2025-09-18T12:43:50.057053Z","shell.execute_reply.started":"2025-09-18T12:43:50.040860Z","shell.execute_reply":"2025-09-18T12:43:50.056383Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"For a better visualisation of the axial CT slice we took as example the whole data set of radiopedia and made a scanning of:\n- Images\n- Masks\n- Overlay (both datasets)\n\nFor a better optimization of resources in this notebook, the .gif's were generated in an outside python enviroment with the following code:\n\n<details>\n  <summary> Click to see the code </summary>\n\n    from pathlib import Path\n    import numpy as np\n    import matplotlib\n    matplotlib.use(\"Agg\")\n    import matplotlib.pyplot as plt\n    import imageio.v2 as imageio\n    \n    # ====== CONFIG ======\n    IMG_PATH = Path(\"images_radiopedia.npy\")\n    MSK_PATH = Path(\"masks_radiopedia.npy\")\n    OUT_DIR  = Path(\"out\")\n    OUT_DIR.mkdir(parents=True, exist_ok=True)\n    \n    # ====== HELPERS ======\n    def load_npy(p: Path) -> np.ndarray:\n        if not p.exists():\n            raise FileNotFoundError(p)\n        arr = np.load(p)\n        return arr\n    \n    def to_hw(img_array: np.ndarray) -> np.ndarray:\n        \"\"\"(N,H,W,1)->(N,H,W), (N,H,W)->(N,H,W).\"\"\"\n        if img_array.ndim == 4 and img_array.shape[-1] == 1:\n            return img_array[..., 0]\n        elif img_array.ndim == 3:\n            return img_array\n        else:\n            raise ValueError(f\"Wait for: (N,H,W,1) o (N,H,W), get: {img_array.shape}\")\n    \n    def masks_to_labels(msk: np.ndarray) -> tuple[np.ndarray, int]:\n        \"\"\"(N,H,W,K)->labels (N,H,W) por argmax; (N,H,W)->labels; return (labels, K).\"\"\"\n        if msk.ndim == 4:\n            labels = np.argmax(msk, axis=-1)\n            K = msk.shape[-1]\n            return labels, K\n        elif msk.ndim == 3:\n            labels = msk\n            K = int(labels.max()) + 1\n            return labels, K\n        else:\n            raise ValueError(f\"Mask unexpected shape: {msk.shape}\")\n    \n    def window_hu(img: np.ndarray, wl=-600, ww=1500) -> np.ndarray:\n        low, high = wl - ww/2, wl + ww/2\n        img_w = np.clip(img, low, high)\n        img_w = (img_w - low) / (high - low)  # 0..1\n        return img_w\n    \n    def make_palette(K: int) -> np.ndarray:\n        base = np.array([\n            [230,  25,  75], [ 60, 180,  75], [  0, 130, 200], [245, 130,  48],\n            [145,  30, 180], [ 70, 240, 240], [240,  50, 230], [210, 245,  60],\n            [250, 190, 190], [  0, 128, 128],\n        ], dtype=np.float32) / 255.0\n        if K <= len(base):\n            return base[:K]\n    \n        rng = np.random.default_rng(42)\n        extra = rng.random((K - len(base), 3))\n        return np.vstack([base, extra])\n    \n    def labels_to_rgb(labels: np.ndarray, palette: np.ndarray) -> np.ndarray:\n        \"\"\"labels (H,W) -> rgb (H,W,3) usando palette[K,3].\"\"\"\n        return palette[labels]\n    \n    def overlay_rgb(base_gray01: np.ndarray, label_map: np.ndarray, palette: np.ndarray, alpha=0.35) -> np.ndarray:\n        \"\"\"base: (H,W) en [0,1], label_map: (H,W) ints, devuelve RGB (H,W,3).\"\"\"\n        base_rgb = np.dstack([base_gray01]*3)\n        color = labels_to_rgb(label_map, palette)\n        # Do not color background class\n        # mask_fg = (label_map != 0)[..., None]\n        mask_fg = np.ones_like(color, dtype=bool)\n        return (1 - alpha)*base_rgb + alpha*color*mask_fg + (~mask_fg)*base_rgb\n    \n    def save_gif(frames: list[np.ndarray], path: Path, fps=15):\n        duration = 1.0 / fps\n        imageio.mimsave(path, [np.clip((f*255).astype(np.uint8), 0, 255) if f.dtype!=np.uint8 else f for f in frames], duration=duration)\n    \n    def save_mp4_matplotlib(frames: list[np.ndarray], path: Path, fps=15):\n        # usa matplotlib + ffmpeg si lo tienes instalado\n        h, w = frames[0].shape[:2]\n        fig = plt.figure(figsize=(w/100, h/100), dpi=100)\n        ax = plt.axes([0,0,1,1]); ax.axis('off')\n        im = ax.imshow(frames[0])\n        import matplotlib.animation as animation\n        ani = animation.ArtistAnimation(fig, [[ax.imshow(fr, animated=True)] for fr in frames], interval=1000//fps, blit=True, repeat=False)\n        ani.save(path, fps=fps, dpi=100)\n        plt.close(fig)\n    \n    # ====== Load data ======\n    imgs_raw = load_npy(IMG_PATH)    # (N,H,W,1) o (N,H,W)\n    imgs = to_hw(imgs_raw)           # -> (N,H,W)\n    N_img = imgs.shape[0]\n    print(\"IMG:\", imgs.shape, imgs.min(), imgs.max(), imgs.mean())\n    \n    labels = None; K = None\n    if MSK_PATH.exists():\n        masks_raw = load_npy(MSK_PATH)    # (N,H,W,K) o (N,H,W)\n        labels, K = masks_to_labels(masks_raw)  # -> (N,H,W), K\n        print(\"MSK labels:\", labels.shape, \"K=\", K)\n    else:\n        print(\"Masks not foundedd: Just will be generated images scanning.\")\n    \n    # ====== (1) Scanning just images ======\n    frames_img = []\n    for i in range(N_img):\n        g = window_hu(imgs[i])        # 0..1\n        rgb = np.dstack([g]*3)        # gris a RGB\n        frames_img.append(rgb)\n    save_gif(frames_img, OUT_DIR/\"scanning_imagen.gif\", fps=15)\n    print(\"Saved:\", OUT_DIR/\"scanning_imagen.gif\")\n    \n    # ====== (2) Scanning just masks ======\n    if labels is not None:\n        N_msk = labels.shape[0]\n        pal = make_palette(K)\n        frames_msk = [labels_to_rgb(labels[i], pal) for i in range(N_msk)]\n        save_gif(frames_msk, OUT_DIR/\"scanning_masks.gif\", fps=15)\n        print(\"Saved:\", OUT_DIR/\"scanning_masks.gif\")\n    \n    # ====== (3) Scanning Overlay (If match) ======\n    if labels is not None:\n        N = min(N_img, labels.shape[0])\n        if N < N_img or N < labels.shape[0]:\n            print(f\"Warning: N missmatch (IMG={N_img}, MSK={labels.shape[0]}). Will be used N={N}.\")\n        pal = make_palette(K)\n        frames_overlay = []\n        for i in range(N):\n            base = window_hu(imgs[i])               # 0..1\n            over = overlay_rgb(base, labels[i], pal, alpha=0.35)\n            frames_overlay.append(np.clip(over, 0, 1))\n        save_gif(frames_overlay, OUT_DIR/\"scanning_overlay.gif\", fps=15)\n        print(\"Saved:\", OUT_DIR/\"barrido_overlay.gif\")\n        \n</details>\n\nAnd next will be showed the overlay .gif","metadata":{}},{"cell_type":"code","source":"from IPython.display import Video, display\n\nvideo_path = \"/kaggle/input/scanning-video/scanning.mp4\"\ndisplay(Video(video_path, embed=True))\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-18T12:43:50.057831Z","iopub.execute_input":"2025-09-18T12:43:50.058181Z","iopub.status.idle":"2025-09-18T12:43:50.627971Z","shell.execute_reply.started":"2025-09-18T12:43:50.058157Z","shell.execute_reply":"2025-09-18T12:43:50.626501Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nCLASS_NAMES = {0: \"GGO\", 1: \"Consolidation\", 2: \"Lungs-other\", 3: \"Background\"}\n\ndef ensure_4ch_masks(m):\n    \"\"\"Return masks as (N,H,W,4) one-hot/binary per channel.\n       Accepts (N,H,W,4) or (N,H,W,2) or (N,H,W) with labels in {0..3}.\n    \"\"\"\n    if m.ndim == 4 and m.shape[-1] == 4:\n        return (m > 0).astype(np.uint8)\n    if m.ndim == 4 and m.shape[-1] == 2:\n        # pad to 4 channels (keep first two, zeros for 2 & 3)\n        m2 = (m > 0).astype(np.uint8)\n        z = np.zeros(m2.shape[:-1] + (2,), dtype=np.uint8)\n        return np.concatenate([m2, z], axis=-1)\n    if m.ndim == 3:\n        # single-channel labels {0..3}\n        N,H,W = m.shape\n        out = np.zeros((N,H,W,4), dtype=np.uint8)\n        for k in range(4):\n            out[..., k] = (m == k).astype(np.uint8)\n        return out\n    raise ValueError(f\"Unexpected mask shape: {m.shape}\")\n\ndef per_split_stats(split_name, masks_4ch):\n    \"\"\"Compute a tidy DataFrame with per-class stats for one split.\"\"\"\n    m = ensure_4ch_masks(masks_4ch)\n    N, H, W, C = m.shape\n    assert C == 4, \"Expecting 4 channels (0..3)\"\n    total_pix = H * W\n\n    # Per-image positives per class\n    pos_img = m.reshape(N, -1, C).sum(axis=1)                 # (N,C)\n    any_img = pos_img > 0                                      # (N,C)\n\n    rows = []\n    for k in range(C):\n        pix_share = m[..., k].mean() * 100.0                   # global pixel share %\n        pct_imgs_with = any_img[:, k].mean() * 100.0           # % images with class\n        pct_empty = (pos_img[:, k] == 0).mean() * 100.0        # % images empty for class\n\n        # Per-image area % for that class\n        area_pct = (pos_img[:, k] / total_pix) * 100.0         # (N,)\n        area_nonzero = area_pct[area_pct > 0]\n        mean_area = float(area_pct.mean()) if N else np.nan\n        median_area = float(np.median(area_nonzero)) if area_nonzero.size else 0.0\n        p90_area = float(np.percentile(area_nonzero, 90)) if area_nonzero.size else 0.0\n\n        rows.append({\n            \"Split\": split_name,\n            \"Class\": CLASS_NAMES.get(k, f\"ch{k}\"),\n            \"Pixel Share % (global)\": round(pix_share, 4),\n            \"% Images with class\": round(pct_imgs_with, 2),\n            \"% Empty masks\": round(pct_empty, 2),\n            \"Mean area % (per image)\": round(mean_area, 4),\n            \"Median area % (non-zero imgs)\": round(median_area, 4),\n            \"P90 area % (non-zero imgs)\": round(p90_area, 4),\n        })\n    return pd.DataFrame(rows)\n\ndef plot_pixel_share(split_name, masks_4ch):\n    m = ensure_4ch_masks(masks_4ch)\n    fracs = m.mean(axis=(0,1,2)) * 100.0\n    labels = [CLASS_NAMES.get(i, f\"ch{i}\") for i in range(m.shape[-1])]\n    plt.figure()\n    plt.bar(range(len(fracs)), fracs)\n    plt.xticks(range(len(fracs)), labels, rotation=10)\n    plt.ylabel(\"Pixel Share (%)\")\n    plt.title(f\"Global pixel share per class — {split_name}\")\n    plt.tight_layout()\n    plt.show()\n\ndef plot_images_with_class(split_name, masks_4ch):\n    m = ensure_4ch_masks(masks_4ch)\n    pos_img = m.reshape(m.shape[0], -1, m.shape[-1]).sum(axis=1)\n    pct = (pos_img > 0).mean(axis=0) * 100.0\n    labels = [CLASS_NAMES.get(i, f\"ch{i}\") for i in range(m.shape[-1])]\n    plt.figure()\n    plt.bar(range(len(pct)), pct)\n    plt.xticks(range(len(pct)), labels, rotation=10)\n    plt.ylabel(\"% of images with class\")\n    plt.title(f\"Presence of classes — {split_name}\")\n    plt.tight_layout()\n    plt.show()\n\ndef plot_area_histograms(split_name, masks_4ch):\n    \"\"\"Histograms of per-image area % for GGO(0) and Cons(1), non-zero images only.\"\"\"\n    m = ensure_4ch_masks(masks_4ch)\n    N, H, W, _ = m.shape\n    total_pix = H * W\n    pos_img = m.reshape(N, -1, 4).sum(axis=1)  # (N,4)\n\n    for k in [0, 1]:  # GGO & Consolidation\n        area = (pos_img[:, k] / total_pix) * 100.0\n        nz = area[area > 0]\n        plt.figure()\n        plt.hist(nz, bins=30)\n        plt.xlabel(\"% of image pixels\")\n        plt.ylabel(\"Count\")\n        plt.title(f\"Lesion area % — {split_name} — {CLASS_NAMES[k]} (non-zero images)\")\n        plt.tight_layout()\n        plt.show()\n\n# ---------- Build the combined stats table ----------\nall_stats = []\n\nif 'masks_medseg' in globals():\n    df_med = per_split_stats(\"MedSeg (train)\", masks_medseg)\n    all_stats.append(df_med)\n    plot_pixel_share(\"MedSeg (train)\", masks_medseg)\n    plot_images_with_class(\"MedSeg (train)\", masks_medseg)\n    plot_area_histograms(\"MedSeg (train)\", masks_medseg)\n\nif 'masks_radiopedia' in globals():\n    df_rad = per_split_stats(\"Radiopaedia (train)\", masks_radiopedia)\n    all_stats.append(df_rad)\n    plot_pixel_share(\"Radiopaedia (train)\", masks_radiopedia)\n    plot_images_with_class(\"Radiopaedia (train)\", masks_radiopedia)\n    plot_area_histograms(\"Radiopaedia (train)\", masks_radiopedia)\n\nif all_stats:\n    stats_df = pd.concat(all_stats, ignore_index=True)\n    display(stats_df)\n    \nelse:\n    print(\"No masks found. Make sure masks_medseg and/or masks_radiopedia are loaded.\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-18T12:43:50.630537Z","iopub.execute_input":"2025-09-18T12:43:50.631011Z","iopub.status.idle":"2025-09-18T12:44:16.977608Z","shell.execute_reply.started":"2025-09-18T12:43:50.630968Z","shell.execute_reply":"2025-09-18T12:44:16.976845Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 2.2 ETL Process and Data Desitions","metadata":{}},{"cell_type":"markdown","source":"With the data analysis and in accordance with the guidelines we determined that for our ETL process:\n\n1) Standardized images\n\n- What: Converted all images to float32, enforced shape (N, 512, 512, 1), and scaled values to [0,1].\n\n- Why: Models train better when inputs share the same scale and shape. This avoids instability from mixed ranges (e.g., 0–255 vs 0–1) and simplifies batching.\n\n2) Prepared masks for the 2 competition classes\n\n- What: From the original 4-channel masks (GGO, Consolidation, Lungs-other, Background), we kept only GGO (0) and Consolidation (1) → masks shaped (N, 512, 512, 2), binary {0,1}.\n\n- Why: The challenge evaluates only these two classes. Dropping the others simplifies the task and metric computation.\n\n3) Prevented “leakage” when splitting Radiopaedia\n\n- What: Radiopaedia contains 9 3D studies (volumes). Without IDs, we inferred volume boundaries by detecting large jumps between consecutive slices. Then we split by complete volumes (some to train, others to validation; ~80/20).\n\n- Why: Slices from the same volume are very similar. Mixing them across train and validation yields overly optimistic validation (the model “recognizes” the study). Volume-wise splitting prevents information leakage.\n\n4) Balanced MedSeg split\n\n- What: For MedSeg, we performed an 80/20 split stratified by has_any_lesion (presence of any GGO/Consolidation). If stratification was impossible (extreme class imbalance), we used a robust fallback that still produces a valid validation set.\n\n- Why: Ensures validation is informative (not nearly all-empty or all-positive) and reflects the test distribution (the test comes from MedSeg).\n\n5) Combined sources for training and validation\n\n- What: Merged the resulting subsets from MedSeg and Radiopaedia to form X_train/Y_train and X_val/Y_val.\n\n- Why: Training with diverse data and validating on a set that includes MedSeg improves generalization toward the official test set.\n\n6) Data augmentation (on-the-fly, not saved)\n\n- What: Defined an on-the-fly augmentation function (flips, small rotation, mild contrast/gamma) applied during training only.\n\n- Why: Shows the model realistic variations, boosting robustness and reducing overfitting—without increasing disk usage (augmented samples are not saved).\n\n7) Saved artifacts\n\nWhy: Makes the workflow reproducible and lets the model start training immediately by loading ready-to-use arrays.","metadata":{}},{"cell_type":"code","source":"# ============================================================\n# COVID-CT Segmentation — Full ETL (saves to /kaggle/working/)\n# ============================================================\n# Steps:\n#  1) Load arrays (if not already loaded) from /kaggle/input/covid-segmentation\n#  2) Select classes (0=GGO, 1=Consolidation) -> masks (N,H,W,2) in {0,1}\n#  3) Normalize images to float32 in [0,1], enforce (N,H,W,1)\n#  4) Infer 9 Radiopaedia volumes (no metadata) via slice-to-slice diffs\n#  5) Split Radiopaedia by VOLUME (~80/20)\n#  6) Split MedSeg with robust stratifier by has_any_lesion (fallbacks if needed)\n#  7) Combine splits (MedSeg + Radiopaedia) into train/val\n#  8) Save to /kaggle/working/prepared\n#  9) Provide on-the-fly augmentation function\n# ============================================================\n\nimport os, json\nimport numpy as np\nimport pandas as pd\nfrom sklearn.model_selection import train_test_split, GroupShuffleSplit\n\n\nOUT_DIR = '/kaggle/working/prepared'\nos.makedirs(OUT_DIR, exist_ok=True)\n\nassert images_radiopedia is not None and masks_radiopedia is not None, \"Radiopaedia arrays not found.\"\nassert images_medseg is not None and masks_medseg is not None, \"MedSeg arrays not found.\"\nassert test_images_medseg is not None, \"MedSeg test array not found.\"\n\n# -----------------\n# Helpers\n# -----------------\ndef ensure_ch1(x: np.ndarray) -> np.ndarray:\n    \"\"\"(N,H,W) or (N,H,W,1) -> float32 (N,H,W,1) in [0,1].\"\"\"\n    x = x.astype(np.float32)\n    if x.ndim == 3:\n        x = x[..., None]\n    xmin, xmax = float(x.min()), float(x.max())\n    if xmin < 0.0 or xmax > 1.0:\n        x = (x - xmin) / (xmax - xmin + 1e-8)\n    return x\n\ndef masks_to_two_classes(m: np.ndarray) -> np.ndarray:\n    \"\"\"Keep only channels 0=GGO and 1=Consolidation -> (N,H,W,2) in {0,1}.\"\"\"\n    if m.ndim == 4:\n        m = (m > 0).astype(np.uint8)\n        if m.shape[-1] >= 2:\n            return m[..., :2]\n        raise ValueError(f\"Mask has {m.shape[-1]} channels; expected >=2.\")\n    if m.ndim == 3:  # label-encoded {0..3}\n        N,H,W = m.shape\n        out = np.zeros((N,H,W,2), dtype=np.uint8)\n        out[...,0] = (m == 0).astype(np.uint8)  # GGO\n        out[...,1] = (m == 1).astype(np.uint8)  # Consolidation\n        return out\n    raise ValueError(f\"Unexpected mask shape: {m.shape}\")\n\ndef has_any_lesion(m2: np.ndarray) -> np.ndarray:\n    \"\"\"Per-slice boolean: any positive pixel in either class.\"\"\"\n    return (m2.reshape(m2.shape[0], -1, 2).sum(axis=1).sum(axis=1) > 0).astype(np.uint8)\n\ndef robust_medseg_split(has_lesion: np.ndarray, test_size=0.2, random_state=42):\n    \"\"\"\n    Returns train_idx, val_idx for MedSeg even if stratification is impossible.\n    Strategy:\n      1) If both classes exist and each has >=2 samples -> stratified split.\n      2) If both classes exist but one has only 1 sample -> force that sample in val; fill rest from majority.\n      3) If only one class exists -> plain random split.\n    \"\"\"\n    N = len(has_lesion)\n    idx = np.arange(N)\n    rng = np.random.default_rng(random_state)\n\n    classes, counts = np.unique(has_lesion, return_counts=True)\n    count_map = dict(zip(classes.tolist(), counts.tolist()))\n    n_pos = count_map.get(1, 0)\n    n_neg = count_map.get(0, 0)\n    n_val = max(1, int(round(test_size * N)))\n\n    # Case 1: stratified possible\n    if n_pos >= 2 and n_neg >= 2:\n        tr_idx, va_idx = train_test_split(\n            idx, test_size=test_size, random_state=random_state, stratify=has_lesion\n        )\n        return tr_idx, va_idx\n\n    pos_idx = idx[has_lesion == 1]\n    neg_idx = idx[has_lesion == 0]\n\n    # Case 2: both classes exist but one has only 1 sample\n    if (n_pos == 1 and n_neg >= 1) or (n_neg == 1 and n_pos >= 1):\n        if n_neg == 1:\n            must_have = neg_idx\n            pool = pos_idx\n        else:\n            must_have = pos_idx\n            pool = neg_idx\n        rest = max(0, n_val - len(must_have))\n        take = rng.choice(pool, size=min(rest, len(pool)), replace=False) if rest > 0 else np.array([], dtype=int)\n        va_idx = np.unique(np.concatenate([must_have, take]))\n        if len(va_idx) < n_val:\n            remaining = np.setdiff1d(idx, va_idx, assume_unique=False)\n            top_up = rng.choice(remaining, size=n_val - len(va_idx), replace=False)\n            va_idx = np.unique(np.concatenate([va_idx, top_up]))\n        tr_idx = np.setdiff1d(idx, va_idx, assume_unique=False)\n        return tr_idx, va_idx\n\n    # Case 3: only one class exists -> random split\n    va_idx = rng.choice(idx, size=n_val, replace=False)\n    tr_idx = np.setdiff1d(idx, va_idx, assume_unique=False)\n    return tr_idx, va_idx\n\ndef infer_radiopaedia_volume_ids(images_: np.ndarray, K: int = 9):\n    \"\"\"\n    Infer K volume blocks from slice-to-slice diffs (no metadata).\n    Returns:\n      volume_id: (N,) int32 in {0..K-1}\n      ranges: list of (start, end) indices per volume\n      diffs: (N-1,) mean abs diffs used for boundary detection\n    \"\"\"\n    X = images_\n    if X.ndim == 4: X = X[...,0]\n    X = X.astype(np.float32)\n    N = len(X)\n    diffs = np.array([np.mean(np.abs(X[i] - X[i-1])) for i in range(1, N)], dtype=np.float32)\n    num_boundaries = K - 1\n    # If N is small or diffs flat, fallback to equal chunks\n    if N <= K or np.allclose(diffs, diffs[0]):\n        chunks = np.array_split(np.arange(N), K)\n        ranges = [(int(c[0]), int(c[-1])+1) for c in chunks]\n    else:\n        boundary_pos = np.argsort(diffs)[-num_boundaries:] + 1\n        boundary_pos = np.sort(boundary_pos)\n        starts = np.concatenate([[0], boundary_pos])\n        ends   = np.concatenate([boundary_pos, [N]])\n        ranges = list(zip(starts, ends))\n        # guard against tiny chunks\n        if any((e - s) < 10 for s, e in ranges):\n            chunks = np.array_split(np.arange(N), K)\n            ranges = [(int(c[0]), int(c[-1])+1) for c in chunks]\n    vol_id = np.empty(N, dtype=np.int32)\n    for vid, (s, e) in enumerate(ranges):\n        vol_id[s:e] = vid\n    return vol_id, ranges, diffs\n\n# -----------------\n# 1) Normalize & select classes\n# -----------------\nX_med = ensure_ch1(images_medseg)               # (100,512,512,1) in [0,1]\nY_med = masks_to_two_classes(masks_medseg)      # (100,512,512,2)\n\nX_rad = ensure_ch1(images_radiopedia)           # (829,512,512,1)\nY_rad = masks_to_two_classes(masks_radiopedia)  # (829,512,512,2)\n\nX_test = ensure_ch1(test_images_medseg)         # (10,512,512,1)\n\n# -----------------\n# 2) Radiopaedia: split by inferred volumes (K=9)\n# -----------------\nvolume_id, vol_ranges, _diffs = infer_radiopaedia_volume_ids(images_radiopedia, K=9)\ngss = GroupShuffleSplit(test_size=0.2, n_splits=1, random_state=42)\n(rad_tr_idx, rad_val_idx), = gss.split(np.arange(len(X_rad)), groups=volume_id)\n\n# -----------------\n# 3) MedSeg: robust stratified split by has_any_lesion (80/20)\n# -----------------\ny_med_has = has_any_lesion(Y_med)\nmed_tr_idx, med_val_idx = robust_medseg_split(y_med_has, test_size=0.2, random_state=42)\n\n# -----------------\n# 4) Combine splits\n# -----------------\nX_train = np.concatenate([X_med[med_tr_idx], X_rad[rad_tr_idx]], axis=0)\nY_train = np.concatenate([Y_med[med_tr_idx], Y_rad[rad_tr_idx]], axis=0)\n\nX_val   = np.concatenate([X_med[med_val_idx], X_rad[rad_val_idx]], axis=0)\nY_val   = np.concatenate([Y_med[med_val_idx], Y_rad[rad_val_idx]], axis=0)\n\nprint(\"Train shapes:\", X_train.shape, Y_train.shape)\nprint(\"Val   shapes:\", X_val.shape,   Y_val.shape)\nval_sources = ([\"MedSeg\"] * len(med_val_idx)) + ([\"Radiopaedia\"] * len(rad_val_idx))\nprint(\"Validation source composition:\\n\", pd.Series(val_sources).value_counts(normalize=True))\n\n# -----------------\n# 5) Save prepared artifacts to /kaggle/working/prepared\n# -----------------\nnp.savez_compressed(os.path.join(OUT_DIR, \"train_arrays.npz\"), X=X_train, Y=Y_train)\nnp.savez_compressed(os.path.join(OUT_DIR, \"val_arrays.npz\"),   X=X_val,   Y=Y_val)\nnp.savez_compressed(os.path.join(OUT_DIR, \"test_medseg.npz\"),  X=X_test)\n\npd.DataFrame({\"rad_train_idx\": rad_tr_idx}).to_csv(os.path.join(OUT_DIR, \"radiopaedia_train_idx.csv\"), index=False)\npd.DataFrame({\"rad_val_idx\":   rad_val_idx}).to_csv(os.path.join(OUT_DIR, \"radiopaedia_val_idx.csv\"),   index=False)\npd.DataFrame({\"med_train_idx\": med_tr_idx}).to_csv(os.path.join(OUT_DIR, \"medseg_train_idx.csv\"),       index=False)\npd.DataFrame({\"med_val_idx\":   med_val_idx}).to_csv(os.path.join(OUT_DIR, \"medseg_val_idx.csv\"),        index=False)\n\nwith open(os.path.join(OUT_DIR, \"info.json\"), \"w\") as f:\n    json.dump({\n        \"images_shape\": list(X_train.shape[1:]),\n        \"masks_shape\":  list(Y_train.shape[1:]),\n        \"train_size\":   int(len(X_train)),\n        \"val_size\":     int(len(X_val)),\n        \"class_map\":    {\"0\": \"GGO\", \"1\": \"Consolidation\"},\n        \"notes\": \"Images normalized to [0,1]; masks are binary (2 channels). Radiopaedia split by inferred volumes (9).\"\n    }, f, indent=2)\n\nprint(f\"Saved artifacts to: {OUT_DIR}\")\nprint(\"Files:\", sorted(os.listdir(OUT_DIR)))\n\n# -----------------\n# 6) On-the-fly augmentations (use inside your training loop)\n# -----------------\nimport cv2\n_rng = np.random.default_rng(123)\n\ndef augment_pair(img: np.ndarray, mask: np.ndarray):\n    \"\"\"\n    img:  (H,W,1) float32 in [0,1]\n    mask: (H,W,2) uint8 {0,1}\n    returns augmented (img, mask)\n    \"\"\"\n    H, W = img.shape[:2]\n\n    # Horizontal flip (p=0.5)\n    if _rng.random() < 0.5:\n        img  = np.flip(img, axis=1).copy()\n        mask = np.flip(mask, axis=1).copy()\n\n    # Vertical flip (p=0.2)\n    if _rng.random() < 0.2:\n        img  = np.flip(img, axis=0).copy()\n        mask = np.flip(mask, axis=0).copy()\n\n    # Small rotation (-12..+12 deg)\n    angle = float(_rng.uniform(-12, 12))\n    M = cv2.getRotationMatrix2D((W/2, H/2), angle, 1.0)\n    img  = cv2.warpAffine(img, M, (W, H), flags=cv2.INTER_LINEAR,  borderMode=cv2.BORDER_REFLECT)\n    m0   = cv2.warpAffine(mask[...,0].astype(np.uint8), M, (W, H), flags=cv2.INTER_NEAREST, borderMode=cv2.BORDER_REFLECT)\n    m1   = cv2.warpAffine(mask[...,1].astype(np.uint8), M, (W, H), flags=cv2.INTER_NEAREST, borderMode=cv2.BORDER_REFLECT)\n    mask = np.stack([m0, m1], axis=-1)\n\n    # Contrast / gamma jitter\n    gamma = float(_rng.uniform(0.85, 1.2))\n    img = np.clip(img, 0, 1) ** gamma\n\n    # Ensure types and dims\n    if img.ndim == 2:\n        img = img[..., None]\n    return img.astype(np.float32), (mask > 0.5).astype(np.uint8)\n\nprint(\"augment_pair(img, mask) ready.\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-18T12:44:16.978585Z","iopub.execute_input":"2025-09-18T12:44:16.978885Z","iopub.status.idle":"2025-09-18T12:45:28.925067Z","shell.execute_reply.started":"2025-09-18T12:44:16.978869Z","shell.execute_reply":"2025-09-18T12:45:28.924447Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\n\ndef has_any_lesion_bin(Y):\n    # Y: (N,H,W,2) uint8 {0,1}\n    return (Y.reshape(Y.shape[0], -1, 2).sum(axis=1).sum(axis=1) > 0)\n\nprint(\"=== SHAPES & DTYPES ===\")\nprint(f\"X_train: {X_train.shape}  dtype={X_train.dtype}  range=({X_train.min():.3f},{X_train.max():.3f})\")\nprint(f\"Y_train: {Y_train.shape}  dtype={Y_train.dtype}  values={{0,1}}\")\nprint(f\"X_val  : {X_val.shape}    dtype={X_val.dtype}    range=({X_val.min():.3f},{X_val.max():.3f})\")\nprint(f\"Y_val  : {Y_val.shape}    dtype={Y_val.dtype}    values={{0,1}}\")\n\nprint(\"\\n=== PER-SET SLICE COUNTS ===\")\nprint(f\"Train slices: {len(X_train)}\")\nprint(f\"Val   slices: {len(X_val)}\")\n\nprint(\"\\n=== LESION PRESENCE (has_any_lesion) ===\")\ntrain_has = has_any_lesion_bin(Y_train)\nval_has   = has_any_lesion_bin(Y_val)\nprint(\"Train: \",\n      pd.Series(train_has.astype(int)).map({0:\"no-lesion\",1:\"lesion\"}).value_counts().to_dict())\nprint(\"Val  : \",\n      pd.Series(val_has.astype(int)).map({0:\"no-lesion\",1:\"lesion\"}).value_counts().to_dict())\n\ntry:\n    train_sources = ([\"MedSeg\"]*len(med_tr_idx)) + ([\"Radiopaedia\"]*len(rad_tr_idx))\n    val_sources   = ([\"MedSeg\"]*len(med_val_idx)) + ([\"Radiopaedia\"]*len(rad_val_idx))\n    print(\"\\n=== SOURCE COMPOSITION ===\")\n    print(\"Train:\", pd.Series(train_sources).value_counts().to_dict())\n    print(\"Val  :\", pd.Series(val_sources).value_counts().to_dict())\nexcept NameError:\n    print(\"\\n(Source indices not available in this scope.)\")\n\nprint(\"\\n=== QUICK SANITY CHECKS ===\")\n# Check few random items for non-empty masks and per-class pixel share\nrng = np.random.default_rng(0)\nfor name, Xs, Ys in [(\"TRAIN\", X_train, Y_train), (\"VAL\", X_val, Y_val)]:\n    idx = rng.choice(len(Xs), size=min(3, len(Xs)), replace=False)\n    for i in idx:\n        ggo_px  = int(Ys[i,...,0].sum())\n        cons_px = int(Ys[i,...,1].sum())\n        total   = Ys[i].shape[0]*Ys[i].shape[1]\n        print(f\"{name} idx={i:4d} | GGO={ggo_px/total*100:.2f}%  Cons={cons_px/total*100:.2f}%  any={ggo_px+cons_px>0}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-18T12:45:28.925812Z","iopub.execute_input":"2025-09-18T12:45:28.926053Z","iopub.status.idle":"2025-09-18T12:45:34.582767Z","shell.execute_reply.started":"2025-09-18T12:45:28.926035Z","shell.execute_reply":"2025-09-18T12:45:34.582000Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os, numpy as np, matplotlib.pyplot as plt\n\ndef _get_arrays():\n    if \"X_train\" in globals() and \"Y_train\" in globals() and \"X_val\" in globals() and \"Y_val\" in globals():\n        return X_train, Y_train, X_val, Y_val\n    prep = \"/kaggle/working/prepared\"\n    tr = np.load(os.path.join(prep, \"train_arrays.npz\"))\n    va = np.load(os.path.join(prep, \"val_arrays.npz\"))\n    return tr[\"X\"], tr[\"Y\"], va[\"X\"], va[\"Y\"]\n\nX_tr, Y_tr, X_va, Y_va = _get_arrays()\n\ndef overlay_masks(img_hw1, mask_hw2, alpha=0.35):\n    \"\"\"\n    img_hw1: (H,W,1) float32 [0,1]\n    mask_hw2: (H,W,2) uint8 {0,1}  [0]=GGO (verde), [1]=Consolidation (rojo)\n    return: (H,W,3) float32 [0,1]\n    \"\"\"\n    img = np.clip(img_hw1[..., 0], 0, 1)\n    H, W = img.shape\n    rgb = np.stack([img, img, img], axis=-1)\n\n    ggo  = mask_hw2[..., 0] > 0\n    cons = mask_hw2[..., 1] > 0\n\n    out = rgb.copy()\n    out[ggo]  = out[ggo]  * (1 - alpha) + np.array([0.0, 1.0, 0.0]) * alpha\n    out[cons] = out[cons] * (1 - alpha) + np.array([1.0, 0.0, 0.0]) * alpha\n    return np.clip(out, 0, 1)\n\ndef show_grid(X, Y, title=\"SET\", n=8, seed=0, cols=4):\n    n = min(n, len(X))\n    idxs = np.random.default_rng(seed).choice(len(X), size=n, replace=False)\n    rows = int(np.ceil(n / cols))\n    plt.figure(figsize=(4*cols, 4*rows))\n    for k, i in enumerate(idxs, 1):\n        ov = overlay_masks(X[i], Y[i], alpha=0.35)\n        ggo_px  = int(Y[i, ..., 0].sum())\n        cons_px = int(Y[i, ..., 1].sum())\n        total   = Y[i].shape[0] * Y[i].shape[1]\n        plt.subplot(rows, cols, k)\n        plt.imshow(ov)\n        plt.title(f\"{title} idx={i} | GGO {100*ggo_px/total:.1f}%  Cons {100*cons_px/total:.1f}%\")\n        plt.axis(\"off\")\n    plt.tight_layout()\n    plt.show()\n\nshow_grid(X_tr, Y_tr, title=\"TRAIN\", n=8, seed=0, cols=4)\nshow_grid(X_va, Y_va, title=\"VAL\",   n=8, seed=1, cols=4)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-18T12:45:34.583605Z","iopub.execute_input":"2025-09-18T12:45:34.583825Z","iopub.status.idle":"2025-09-18T12:45:36.878427Z","shell.execute_reply.started":"2025-09-18T12:45:34.583799Z","shell.execute_reply":"2025-09-18T12:45:36.877503Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"During ETL we never change the spatial grid or the pairing logic—only the mask channels:\n\n- **Same index, same slice:** We construct X_* (images) and Y_* (masks) using the same slice indices for each split/source. Thus X_train[i] and Y_train[i] (and likewise for validation) come from the exact same CT slice.\n\n- **Same resolution and geometry:** Both images and masks keep the original 512×512 grid. We do not resample, crop, or pad in the ETL; the only operation on masks is channel selection (keep 0=GGO, 1=Consolidation; drop the other two). This preserves all pixel coordinates.\n\n- **No independent shuffling:** We split by indices (and, for Radiopaedia, by inferred volume groups) and then apply those indices simultaneously to images and masks. There is never a reordering of images without applying the same reordering to masks.\n\n- **Identical preprocessing path:** Images are normalized to [0,1] and coerced to (H,W,1); masks are cast to binary (H,W,2) by channel slicing. No geometric transforms occur in ETL, so there is no chance to misalign image and mask.","metadata":{}},{"cell_type":"markdown","source":"## 3. Model","metadata":{}},{"cell_type":"markdown","source":"### 3.1 Why U-Net for this task\n\nAfter consider between 3 sementation models (U-Net, DeepLab and SegNet), we decide use U-Net, desgined in 2015 by Ronnenber O., Fischer P. and Brox T. at the University of Freiburg, Germany. It follows a symmetric encoder–decoder (“U-shaped”) design with skip connections that copy high-resolution features from the encoder to the decoder, enabling precise boundary localization even with limited training data.U-Net were used in order to resolve a segmentation problem, similar to our competition, U-net result with some advatages like:\n\n- **Built for medical imaging:** U-Net was designed for biomedical segmentation and is known to work well with small labeled datasets by leveraging strong skip connections and data augmentation. \n\n- **Simple, fast, reliable:** It’s easy to train end-to-end, delivers sharp boundaries via encoder–decoder skips, and has a very large, supportive community/ecosystem. \n\n- **Good fit for our data scale:** Compared with heavier designs (e.g., DeepLab), plain U-Net typically attains strong performance on limited medical data with modest compute.\n\n\n### 3.2 Our architecture configuration\n\n1) Input & preprocessing\n\n    - What: CT input (1, H, W) scaled to [0,1]; train at 256×256 (from ETL). Validation/inference can be 512×512 if the model/VRAM allow.\n    \n    - Why: 256×256 speeds training and fits VRAM; keep a consistent pipeline across train/val/test.\n    \n    - Alternatives: Train directly at 512×512 (more detail, smaller batch) or use patching.\n\n2) Depth (4-down / 4-up)\n\n    - What: 4 encoder stages with 2×2 max-pool (down to 1/16 res) and 4 decoder stages with skip connections.\n    \n    - Why: Good trade-off between context and detail for 256×256 / 512×512 with moderate VRAM.\n    \n    - Alternatives: 3-down (lighter, less context) or 5-down (more capacity/VRAM).\n\n3) Channel width (64 → 1024 → 64)\n\n    - What: Stage widths: 64, 128, 256, 512; bottleneck 1024; mirrored in decoder.\n    \n    - Why: The classic and stable recipe for small medical datasets.\n    \n    - Alternatives: Halved widths (32–512) to save VRAM; larger widths if data/compute allow.\n\n4) Basic block: Conv3×3 + BatchNorm + ReLU (×2 per stage)\n\n    - What: Two 3×3, stride=1, padding=1 convolutions with BN and ReLU at each level.\n    \n    - Why: Proven pattern—stable training and sharp boundaries.\n    \n    - Alternatives: Instance/GroupNorm (very small batches), LeakyReLU/SiLU, residual (ResU-Net).\n\n5) Downsampling: MaxPool 2×2 (stride 2)\n\n    - What: Resolution reduction after each encoder block (except the bottom).\n    \n    - Why: Simple, robust; preserves scale semantics.\n    \n    - Alternatives: Strided conv (s=2) or AvgPool.\n\n6) Upsampling: Bilinear ×2 + Conv2×2, with concat skip\n\n    - What: ×2 interpolation, concatenate encoder skip, then two Conv+BN+ReLU.\n    \n    - Why: Avoids checkerboard artifacts, lightweight, works very well.\n    \n    - Alternatives: TransposedConv 2×2 (learned upsampling; slightly heavier).\n\n7) Output head: Conv1×1 → 2 logits, Sigmoid at inference\n\n    - What: Map 64 channels to 2 (GGO, Consolidation); apply per-channel Sigmoid at inference.\n    \n    - Why: This is a per-pixel multilabel problem (classes are not mutually exclusive).\n    \n    - Alternatives: Softmax for mutually exclusive classes (not our case).\n\n8) Loss function: BCEWithLogits + Dice (mean over channels)\n\n    - What: Combine BCE (on logits) with Dice.\n    \n    - Why: Dice copes with class imbalance; BCE stabilizes gradients.\n    \n    - Alternatives: Tversky (α,β) to bias FN/FP, Focal (often with Dice) for many negatives.\n\n9) Optimizer & learning rate\n\n    - What: AdamW (lr 1e-3, weight decay 1e-4), cosine annealing LR (T_max = epochs).\n    \n    - Why: Fast convergence and good generalization; cosine is simple and effective.\n    \n    - Alternatives: Adam (no WD), ReduceLROnPlateau (patience 5, factor 0.5), SGD+momentum (better with lots of data).\n  \n### 3.3 Regulatizarion Techniques\n\nOur implementation of the model, we applied L2 regularization, also known as ridge regression through the optimizer AdamW from PyTorch. Code line:\noptimizer = torch.optim.AdamW(model.parameters(), lr=1e-3, weight_decay=1e-4)\nWithin those parameters, we add the parameters lr = 1e-3 (learning rate) and weight_decay=1e-4.  This technique penalizes large weight values by calculating the square of the weights during training, encouraging the network to maintain smaller and more stable parameters, and reducing the risk of overfitting.\nIn our implementation, the learning_rate parameter is established to control the size of the weight training updates. According to the PyTorch AdamW documentation, lr = 1e-3 is the default setting that works best for the optimizer, so we chose not to move it. As a higher learning rate speed would speed up the training but make it unstable for our model. And a lower training rate would make the model more stable, but slower at learning.\n\nwe chose our parameter 1e-4 for weight_decay, which can be tuned to control the strength coefficient of the weight penalties. A smaller value (1e-5) would cause weaker regularization and not show much change, while a larger value (1e-3), would apply stronger penalties that would negatively affect the model. We had to tune the parameter closely to ensure that it reduced overfitting while not being too high. This balance contributed to stabilizing the training process and improving validation performance in our U-Net model.\n\nAlso we use augmentation as showed in the previous blocks, this is part of a regularization technique because it helps to prevent overfitting increasing the divertsity of data by an artifitial samples generation, changing charactetistics of the data like rotatiom, changing color, contrast, illumination or size characteristics.","metadata":{}},{"cell_type":"markdown","source":"## 4. Training\n\n### 4.1 Hardware\n\nOur model were trained with the buil-in Kaggle resources:\n- CPU: Intel Xeon 2.20 GHz, 4 vCPU cores.\n- RAM: 32 GB.\n- GPU: NVIDIA Tesla P100 GPU, 3584 Cuda cores, 16 GB.\n\nWe the available resources we get a run time of 2 hour approximately.\n\n### 4.2 U-Net configuration\n\n\n- SEED: 42\n\n- Expected shapes: img: (H,W,1), mask: (H,W,2)\n\n- Augmentations: augment_pair method (declared above).\n\n- Batch size: 8\n\n- num_workers: 2\n\n- Input channels: 1\n\n- Classes: 2\n\n- Depth: classic 4 downs / 4 ups\n\n- Down path: 64 → 128 → 256 → 512 → bottleneck 1024\n\n- Up path with skip-concat: (1024+512)→512, (512+256)→256, (256+128)→128, (128+64)→64\n\n- Upsampling: bilinear upsample (align_corners=False)\n\n- Blocks: Conv(3×3, bias=False) + BatchNorm2d + ReLU ×2 per block\n\n- Training loss: BCEWithLogitsLoss + soft Dice (per-class, eps 1e-7)\n\n- Optimizer: AdamW(lr=1e-3, weight_decay=1e-4)\n(betas default to PyTorch: (0.9, 0.999), eps=1e-8)\n\n- Epochs: 60\n\n- Early stop: not enabled; trains full 60 epochs while tracking best","metadata":{}},{"cell_type":"code","source":"#%%skip_if True\n# ============================\n# U-Net (2D) + Training & Eval\n# ============================\nimport os, math, time, random, contextlib\nimport numpy as np\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nfrom torch.utils.data import Dataset, DataLoader\n\n# --------------------\n# Reproducibility\n# --------------------\nSEED = 42\nrandom.seed(SEED); np.random.seed(SEED); torch.manual_seed(SEED)\ntorch.backends.cudnn.deterministic = True\ntorch.backends.cudnn.benchmark = False\n\n# --------------------\n# Device & AMP (new API; only if CUDA)\n# --------------------\ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\nUSE_CUDA = (device.type == \"cuda\")\nprint(\"Device:\", device)\n\nif USE_CUDA:\n    scaler = torch.amp.GradScaler('cuda')\n    autocast_ctx = lambda: torch.amp.autocast('cuda', dtype=torch.float16)\nelse:\n    scaler = None\n    autocast_ctx = contextlib.nullcontext  # no-op on CPU\n\n# --------------------\n# Dataset (uses your augment_pair from ETL if provided)\n# --------------------\nclass NumpySegDataset(Dataset):\n    def __init__(self, X, Y=None, train=False, apply_aug_fn=None):\n        self.X = X; self.Y = Y\n        self.train = train; self.apply_aug_fn = apply_aug_fn\n    def __len__(self): return len(self.X)\n    def __getitem__(self, i):\n        img = self.X[i]  # (H,W,1)\n        if self.Y is None:\n            if self.apply_aug_fn is not None and self.train:\n                img, _ = self.apply_aug_fn(img, np.zeros((*img.shape[:2], 2), dtype=np.uint8))\n            x = torch.from_numpy(img).permute(2,0,1).float()\n            return x\n        mask = self.Y[i] # (H,W,2)\n        if self.apply_aug_fn is not None and self.train:\n            img, mask = self.apply_aug_fn(img, mask)\n        x = torch.from_numpy(img).permute(2,0,1).float()\n        y = torch.from_numpy(mask).permute(2,0,1).float()\n        return x, y\n\n# --------------------\n# U-Net (classic, 4-down/4-up, 64..1024)\n# --------------------\nclass DoubleConv(nn.Module):\n    def __init__(self, in_ch, out_ch):\n        super().__init__()\n        self.block = nn.Sequential(\n            nn.Conv2d(in_ch, out_ch, 3, padding=1, bias=False),\n            nn.BatchNorm2d(out_ch),\n            nn.ReLU(inplace=True),\n            nn.Conv2d(out_ch, out_ch, 3, padding=1, bias=False),\n            nn.BatchNorm2d(out_ch),\n            nn.ReLU(inplace=True),\n        )\n    def forward(self, x): return self.block(x)\n\nclass Down(nn.Module):\n    def __init__(self, in_ch, out_ch):\n        super().__init__()\n        self.pool = nn.MaxPool2d(2)\n        self.conv = DoubleConv(in_ch, out_ch)\n    def forward(self, x): return self.conv(self.pool(x))\n\nclass Up(nn.Module):\n    def __init__(self, in_ch, out_ch, bilinear=True):\n        super().__init__()\n        if bilinear:\n            self.up = nn.Upsample(scale_factor=2, mode='bilinear', align_corners=False)\n            self.conv = DoubleConv(in_ch, out_ch)\n        else:\n            self.up = nn.ConvTranspose2d(in_ch//2, in_ch//2, 2, stride=2)\n            self.conv = DoubleConv(in_ch, out_ch)\n    def forward(self, x1, x2):\n        x1 = self.up(x1)\n        dy = x2.size(2) - x1.size(2)\n        dx = x2.size(3) - x1.size(3)\n        x1 = F.pad(x1, [dx//2, dx - dx//2, dy//2, dy - dy//2])\n        x = torch.cat([x2, x1], dim=1)\n        return self.conv(x)\n\nclass OutConv(nn.Module):\n    def __init__(self, in_ch, out_ch):\n        super().__init__()\n        self.conv = nn.Conv2d(in_ch, out_ch, kernel_size=1)\n    def forward(self, x): return self.conv(x)\n\nclass UNet2D(nn.Module):\n    def __init__(self, n_channels=1, n_classes=2, bilinear=True):\n        super().__init__()\n        self.inc   = DoubleConv(n_channels, 64)\n        self.down1 = Down(64, 128)\n        self.down2 = Down(128, 256)\n        self.down3 = Down(256, 512)\n        self.down4 = Down(512, 1024)\n        self.up1   = Up(1024+512, 512, bilinear)\n        self.up2   = Up(512+256, 256, bilinear)\n        self.up3   = Up(256+128, 128, bilinear)\n        self.up4   = Up(128+64,  64,  bilinear)\n        self.outc  = OutConv(64, n_classes)\n        for m in self.modules():\n            if isinstance(m, nn.Conv2d):\n                nn.init.kaiming_normal_(m.weight)\n            elif isinstance(m, (nn.BatchNorm2d, nn.GroupNorm)):\n                nn.init.ones_(m.weight); nn.init.zeros_(m.bias)\n    def forward(self, x):\n        x1 = self.inc(x)\n        x2 = self.down1(x1)\n        x3 = self.down2(x2)\n        x4 = self.down3(x3)\n        x5 = self.down4(x4)\n        x  = self.up1(x5, x4)\n        x  = self.up2(x,  x3)\n        x  = self.up3(x,  x2)\n        x  = self.up4(x,  x1)\n        return self.outc(x)  # (B,2,H,W)\n\n# --------------------\n# Loss: BCEWithLogits + Dice\n# --------------------\nclass BCEDiceLoss(nn.Module):\n    def __init__(self, eps=1e-7):\n        super().__init__()\n        self.bce = nn.BCEWithLogitsLoss()\n        self.eps = eps\n    def dice_soft(self, logits, targets):\n        probs = torch.sigmoid(logits)\n        dims = (0,2,3)\n        inter = (probs * targets).sum(dims)\n        den   = (probs + targets).sum(dims).clamp_min(self.eps)\n        dice  = (2*inter) / den\n        return 1 - dice.mean()\n    def forward(self, logits, targets):\n        return self.bce(logits, targets) + self.dice_soft(logits, targets)\n\n# --------------------\n# DataLoaders (hooks your arrays & augment_pair)\n# --------------------\ntry:\n    X_train, Y_train, X_val, Y_val, X_test  # noqa\nexcept NameError:\n    prep = \"/kaggle/working/prepared\"\n    tr = np.load(os.path.join(prep, \"train_arrays.npz\"))\n    va = np.load(os.path.join(prep, \"val_arrays.npz\"))\n    te = np.load(os.path.join(prep, \"test_medseg.npz\"))\n    X_train, Y_train = tr[\"X\"], tr[\"Y\"]\n    X_val,   Y_val   = va[\"X\"], va[\"Y\"]\n    X_test          = te[\"X\"]\n    print(\"Loaded arrays from prepared/.\")\n\nif \"augment_pair\" not in globals():\n    def augment_pair(img, mask): return img, mask  # identity\n\nBATCH = 8\ntrain_ds = NumpySegDataset(X_train, Y_train, train=True,  apply_aug_fn=augment_pair)\nval_ds   = NumpySegDataset(X_val,   Y_val,   train=False, apply_aug_fn=None)\npin = USE_CUDA  # only pin on CUDA\ntrain_dl = DataLoader(train_ds, batch_size=BATCH, shuffle=True,  num_workers=2, pin_memory=pin)\nval_dl   = DataLoader(val_ds,   batch_size=BATCH, shuffle=False, num_workers=2, pin_memory=pin)\n\n# --------------------\n# Model, optimizer, schedule\n# --------------------\nmodel = UNet2D(n_channels=1, n_classes=2, bilinear=True).to(device)\ncriterion = BCEDiceLoss()\noptimizer = torch.optim.AdamW(model.parameters(), lr=1e-3, weight_decay=1e-4)\nEPOCHS = 60\nscheduler = torch.optim.lr_scheduler.CosineAnnealingLR(optimizer, T_max=EPOCHS)\n\n# --------------------\n# Metrics\n# --------------------\n@torch.no_grad()\ndef dice_score(logits, targets, thr=0.5, eps=1e-7):\n    probs = torch.sigmoid(logits)\n    preds = (probs > thr).float()\n    dims = (0,2,3)\n    inter = (preds * targets).sum(dims)\n    den   = (preds + targets).sum(dims).clamp_min(eps)\n    dice_per_ch = (2*inter) / den\n    return dice_per_ch.detach().cpu().numpy(), float(dice_per_ch.mean().item())\n\n# =====================\n# Metrics + Viz helpers\n# =====================\nimport copy, matplotlib.pyplot as plt, pandas as pd\n\nTHR = 0.45  # threshold for binarizing probs (tune later)\n\n@torch.no_grad()\ndef evaluate_epoch(model, loader, criterion, thr=0.5, device=device, use_cuda=USE_CUDA):\n    \"\"\"\n    Returns:\n      loss_avg, dice_mean, acc_mean, conf (dict with 2x2 matrices per class)\n    Conf matrices are in sklearn order: [[TN, FP],[FN, TP]]\n    \"\"\"\n    model.eval()\n    n_pix = 0\n    loss_sum = 0.0\n\n    # Aggregates for Dice\n    inter_sum = np.zeros(2, dtype=np.float64)\n    den_sum   = np.zeros(2, dtype=np.float64)\n\n    # Aggregates for accuracy & confusion\n    TP = np.zeros(2, dtype=np.float64)\n    FP = np.zeros(2, dtype=np.float64)\n    FN = np.zeros(2, dtype=np.float64)\n    TN = np.zeros(2, dtype=np.float64)\n\n    for xb, yb in loader:\n        xb = xb.to(device, non_blocking=use_cuda)\n        yb = yb.to(device, non_blocking=use_cuda)  # (B,2,H,W)\n\n        logits = model(xb)\n        loss = criterion(logits, yb)\n        loss_sum += loss.item() * xb.size(0)\n\n        probs = torch.sigmoid(logits)\n        preds = (probs > thr).float()\n\n        # Dice aggregations\n        inter = (preds * yb).sum(dim=(0,2,3))       # (2,)\n        den   = (preds + yb).sum(dim=(0,2,3))       # (2,)\n        inter_sum += inter.detach().cpu().numpy()\n        den_sum   += den.detach().cpu().numpy()\n\n        # Accuracy + confusion per class (binary, per pixel)\n        for c in range(2):\n            p = preds[:, c]\n            t = yb[:, c]\n            TP[c] += (p * t).sum().item()\n            FP[c] += (p * (1 - t)).sum().item()\n            FN[c] += ((1 - p) * t).sum().item()\n            TN[c] += ((1 - p) * (1 - t)).sum().item()\n\n        n_pix += yb.numel() // 2  # total pixels per class summed over batch\n\n    # Dice per-class and mean\n    dice_pc = (2 * inter_sum + 1e-7) / (den_sum + 1e-7)\n    dice_mean = float(np.nanmean(dice_pc))\n\n    # Pixel accuracy per class and mean\n    acc_pc = (TP + TN) / (TP + FP + FN + TN + 1e-7)\n    acc_mean = float(np.nanmean(acc_pc))\n\n    loss_avg = loss_sum / len(loader.dataset)\n\n    conf = {\n        \"GGO\": np.array([[TN[0], FP[0]], [FN[0], TP[0]]], dtype=np.float64),\n        \"CONS\": np.array([[TN[1], FP[1]], [FN[1], TP[1]]], dtype=np.float64),\n        \"dice_pc\": dice_pc, \"acc_pc\": acc_pc,\n    }\n    return loss_avg, dice_mean, acc_mean, conf\n\ndef plot_training_curves(history, out_dir=\"/kaggle/working\"):\n    epochs = history[\"epoch\"]\n    # Accuracy\n    plt.figure(figsize=(7,5))\n    plt.plot(epochs, history[\"train_acc\"], label=\"Train acc\")\n    plt.plot(epochs, history[\"val_acc\"],   label=\"Val acc\")\n    plt.xlabel(\"Epoch\"); plt.ylabel(\"Pixel accuracy (mean over classes)\")\n    plt.title(\"Accuracy — train vs val\")\n    plt.legend(); plt.tight_layout()\n    plt.savefig(f\"{out_dir}/accuracy.png\", dpi=150)\n    plt.show()\n\n    # Loss\n    plt.figure(figsize=(7,5))\n    plt.plot(epochs, history[\"train_loss\"], label=\"Train loss\")\n    plt.plot(epochs, history[\"val_loss\"],   label=\"Val loss\")\n    plt.xlabel(\"Epoch\"); plt.ylabel(\"Loss\")\n    plt.title(\"Loss — train vs val\")\n    plt.legend(); plt.tight_layout()\n    plt.savefig(f\"{out_dir}/loss.png\", dpi=150)\n    plt.show()\n\ndef plot_confusions(conf, out_path=\"/kaggle/working/confusion.png\"):\n    fig, axs = plt.subplots(1, 2, figsize=(10,4))\n    titles = [\"GGO\", \"Consolidation\"]\n    mats = [conf[\"GGO\"], conf[\"CONS\"]]\n    for ax, title, M in zip(axs, titles, mats):\n        im = ax.imshow(M, interpolation=\"nearest\")\n        ax.set_title(f\"{title} — Confusion\")\n        ax.set_xticks([0,1]); ax.set_yticks([0,1])\n        ax.set_xticklabels([\"TN\",\"FP\"]); ax.set_yticklabels([\"FN\",\"TP\"])\n        # annotate\n        for i in range(2):\n            for j in range(2):\n                ax.text(j, i, f\"{int(M[i,j])}\", ha=\"center\", va=\"center\")\n    fig.tight_layout()\n    plt.savefig(out_path, dpi=150)\n    plt.show()\n\n# ======================================\n# Train / Validate\n# ======================================\nbest_val = -1.0\nckpt_path = \"/kaggle/working/best_unet.pt\"\n\n# extra loader to evaluate TRAIN without augmentations\ntrain_eval_ds = NumpySegDataset(X_train, Y_train, train=False, apply_aug_fn=None)\ntrain_eval_dl = DataLoader(train_eval_ds, batch_size=BATCH, shuffle=False,\n                           num_workers=2, pin_memory=USE_CUDA)\n\nhistory = {\"epoch\": [], \"train_loss\": [], \"val_loss\": [],\n           \"train_dice\": [], \"val_dice\": [],\n           \"train_acc\": [],  \"val_acc\": []}\n\nfor epoch in range(1, EPOCHS+1):\n    # ------- TRAIN -------\n    model.train()\n    tr_loss = 0.0\n    for xb, yb in train_dl:\n        xb = xb.to(device, non_blocking=USE_CUDA)\n        yb = yb.to(device, non_blocking=USE_CUDA)\n        optimizer.zero_grad(set_to_none=True)\n        with autocast_ctx():\n            logits = model(xb)\n            loss = criterion(logits, yb)\n        if USE_CUDA:\n            scaler.scale(loss).backward()\n            scaler.step(optimizer)\n            scaler.update()\n        else:\n            loss.backward()\n            optimizer.step()\n        tr_loss += loss.item() * xb.size(0)\n    tr_loss /= len(train_ds)\n\n    # ------- EVAL (train & val) -------\n    train_loss_eval, train_dice, train_acc, _      = evaluate_epoch(model, train_eval_dl, criterion, thr=THR)\n    val_loss_eval,   val_dice,   val_acc,   val_c  = evaluate_epoch(model, val_dl,       criterion, thr=THR)\n\n    # step scheduler (keep your cosine)\n    scheduler.step()\n\n    # track best on val Dice\n    if val_dice > best_val:\n        best_val = val_dice\n        torch.save(model.state_dict(), ckpt_path)\n\n    # ---- record & print ----\n    history[\"epoch\"].append(epoch)\n    history[\"train_loss\"].append(tr_loss)\n    history[\"val_loss\"].append(val_loss_eval)\n    history[\"train_dice\"].append(train_dice)\n    history[\"val_dice\"].append(val_dice)\n    history[\"train_acc\"].append(train_acc)\n    history[\"val_acc\"].append(val_acc)\n\n    print(f\"Epoch {epoch:03d}/{EPOCHS} | lr={scheduler.get_last_lr()[0]:.2e} | \"\n          f\"train_loss={tr_loss:.4f} val_loss={val_loss_eval:.4f} \"\n          f\"train_acc={train_acc:.4f} val_acc={val_acc:.4f} \"\n          f\"train_dice={train_dice:.4f} val_dice={val_dice:.4f} best={best_val:.4f}\")\n\n# Save final history & plots\nhist_df = pd.DataFrame(history)\nhist_csv = \"/kaggle/working/history.csv\"\nhist_df.to_csv(hist_csv, index=False)\nprint(\"Saved training history to:\", hist_csv)\nplot_training_curves(history)\nbest_state = torch.load(ckpt_path, map_location=device)\nmodel.load_state_dict(best_state)\n_, _, _, val_conf = evaluate_epoch(model, val_dl, criterion, thr=THR)\nplot_confusions(val_conf, out_path=\"/kaggle/working/confusion.png\")\n\n\nprint(\"Best val Dice:\", best_val)\nprint(\"Saved best weights to:\", ckpt_path)\nprint(\"Figures saved to: /kaggle/working/accuracy.png, /kaggle/working/loss.png, /kaggle/working/confusion.png\")","metadata":{"execution":{"iopub.status.busy":"2025-09-18T12:45:36.879361Z","iopub.execute_input":"2025-09-18T12:45:36.879590Z","iopub.status.idle":"2025-09-18T12:45:36.889269Z","shell.execute_reply.started":"2025-09-18T12:45:36.879571Z","shell.execute_reply":"2025-09-18T12:45:36.888502Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 5. Results","metadata":{}},{"cell_type":"code","source":"#%%skip_if True\n\nimport numpy as np\nimport torch\n\n@torch.no_grad()\ndef eval_metrics(val_loader, thr=0.45, eps=1e-7, count_empty_as=1.0):\n    \n    model.eval()\n\n    inter = np.zeros(2, dtype=np.float64)\n    union = np.zeros(2, dtype=np.float64)\n    den_dice = np.zeros(2, dtype=np.float64)\n    tp = np.zeros(2, dtype=np.float64)\n    fp = np.zeros(2, dtype=np.float64)\n    fn = np.zeros(2, dtype=np.float64)\n    tn = np.zeros(2, dtype=np.float64)\n\n    empty_cases = np.zeros(2, dtype=np.int64)\n    used_cases  = np.zeros(2, dtype=np.int64)\n\n    for xb, yb in val_loader:\n        xb = xb.to(device); yb = yb.to(device)           # yb: (B,2,H,W) {0,1}\n        logits = model(xb)\n        probs = torch.sigmoid(logits)\n        preds = (probs > thr).float()\n\n        for c in range(2):\n            p = preds[:, c]\n            t = yb[:, c]\n\n            tp[c] += (p * t).sum().item()\n            fp[c] += (p * (1 - t)).sum().item()\n            fn[c] += ((1 - p) * t).sum().item()\n            tn[c] += ((1 - p) * (1 - t)).sum().item()\n\n            inter_c = (p * t).sum(dim=(1,2))\n            sum_c   = (p + t).sum(dim=(1,2))\n            union_c = ((p + t) > 0).sum(dim=(1,2))\n            both_empty = (sum_c == 0)\n\n            if count_empty_as is None:\n                mask = ~both_empty\n                inter[c]   += inter_c[mask].sum().item()\n                den_dice[c] += (sum_c[mask]).sum().item()\n                union[c]   += ((p[mask] + t[mask]) > 0).sum().item()\n                used_cases[c] += int(mask.sum().item())\n            else:\n                inter[c]   += inter_c.sum().item()\n                den_dice[c] += sum_c.sum().item()\n                union[c]   += ((p + t) > 0).sum().item()\n                empty_cases[c] += int(both_empty.sum().item())\n                used_cases[c]  += int(len(both_empty))\n\n    # Metrics per class\n    dice = np.zeros(2, dtype=np.float64)\n    iou  = np.zeros(2, dtype=np.float64)\n    acc  = np.zeros(2, dtype=np.float64)\n\n    for c in range(2):\n        if count_empty_as is None:\n            dice[c] = (2*inter[c] + eps) / (den_dice[c] + eps) if used_cases[c] > 0 else np.nan\n        else:\n            dice[c] = (2*inter[c] + eps) / (den_dice[c] + eps)\n\n        iou[c] = dice[c] / (2 - dice[c]) if np.isfinite(dice[c]) else np.nan\n\n        # Accuracy per píxel\n        total = tp[c] + fp[c] + fn[c] + tn[c]\n        acc[c] = (tp[c] + tn[c]) / total if total > 0 else np.nan\n\n    metrics = {\n        \"dice_ggo\": float(dice[0]), \"dice_cons\": float(dice[1]),\n        \"iou_ggo\": float(iou[0]),   \"iou_cons\":  float(iou[1]),\n        \"acc_ggo\": float(acc[0]),   \"acc_cons\":  float(acc[1]),\n        \"dice_mean\": float(np.nanmean(dice)),\n        \"iou_mean\":  float(np.nanmean(iou)),\n        \"acc_mean\":  float(np.nanmean(acc)),\n        \"empty_cases_per_class\": {\"ggo\": int(empty_cases[0]), \"cons\": int(empty_cases[1])},\n        \"used_images_per_class\": {\"ggo\": int(used_cases[0]), \"cons\": int(used_cases[1])},\n        \"thr\": float(thr),\n        \"confusion\": {\n            \"ggo\": {\"tp\": int(tp[0]), \"fp\": int(fp[0]), \"fn\": int(fn[0]), \"tn\": int(tn[0])},\n            \"cons\":{\"tp\": int(tp[1]), \"fp\": int(fp[1]), \"fn\": int(fn[1]), \"tn\": int(tn[1])},\n        }\n    }\n    return metrics\n\ndef print_metrics_table(metrics, save_csv=\"/kaggle/working/val_metrics_summary.csv\"):\n    # Filas per class\n    rows = []\n    rows.append({\n        \"class\": \"ggo\",\n        \"dice\": metrics[\"dice_ggo\"],\n        \"iou\":  metrics[\"iou_ggo\"],\n        \"acc\":  metrics[\"acc_ggo\"],\n    })\n    rows.append({\n        \"class\": \"cons\",\n        \"dice\": metrics[\"dice_cons\"],\n        \"iou\":  metrics[\"iou_cons\"],\n        \"acc\":  metrics[\"acc_cons\"],\n    })\n    # Mean Rows\n    rows.append({\n        \"class\": \"mean\",\n        \"dice\": metrics[\"dice_mean\"],\n        \"iou\":  metrics[\"iou_mean\"],\n        \"acc\":  metrics[\"acc_mean\"],\n    })\n\n    df = pd.DataFrame(rows, columns=[\"class\",\"dice\",\"iou\",\"acc\"])\n\n    df_fmt = df.copy()\n    for col in [\"dice\",\"iou\",\"acc\"]:\n        df_fmt[col] = df_fmt[col].apply(lambda x: np.nan if x is None else float(x))\n        df_fmt[col] = df_fmt[col].map(lambda v: f\"{v:.3f}\" if pd.notnull(v) else \"nan\")\n\n    thr = metrics.get(\"thr\", None)\n    print(\"=\"*46)\n    print(f\" Validation metrics table (thr={thr}) \".center(46, \"=\"))\n    print(\"=\"*46)\n    print(df_fmt.to_string(index=False))\n    print(\"=\"*46)\n\n    df.to_csv(save_csv, index=False)\n    print(\"Saved table CSV ->\", save_csv)\n\nmetrics = eval_metrics(val_dl, thr=0.45, count_empty_as=1.0)\nprint_metrics_table(metrics)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-18T12:45:36.890096Z","iopub.execute_input":"2025-09-18T12:45:36.890316Z","iopub.status.idle":"2025-09-18T12:45:36.907010Z","shell.execute_reply.started":"2025-09-18T12:45:36.890302Z","shell.execute_reply":"2025-09-18T12:45:36.906237Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#%%skip_if True\n\nimport os, cv2, matplotlib.pyplot as plt\nimport numpy as np\nimport torch\n\nbest_state = torch.load(ckpt_path, map_location=device)\nmodel.load_state_dict(best_state)\nmodel.eval()\n\n\ndef render_examples(val_loader, out_dir=\"/kaggle/working/figs\", n_samples=8, thr=0.45, show=True):\n    os.makedirs(out_dir, exist_ok=True)\n    model.eval()\n    saved = 0\n\n    with torch.no_grad():\n        for xb, yb in val_loader:\n            xb = xb.to(device); yb = yb.to(device)\n            probs = torch.sigmoid(model(xb))\n            preds = (probs > thr).float()\n\n            for i in range(xb.size(0)):\n                if saved >= n_samples:\n                    return\n                img = xb[i,0].cpu().numpy()                       # (H,W)\n                gt  = yb[i].cpu().numpy().astype(np.uint8)        # (2,H,W)\n                pr  = preds[i].cpu().numpy().astype(np.uint8)     # (2,H,W)\n\n                # Dice por imagen/clase\n                def dice_np(p, t, eps=1e-7):\n                    inter = (p & t).sum()\n                    den   = (p | t).sum() if (p.max()==1 and t.max()==1) else (p + t).sum()\n                    return (2*inter + eps)/(den + eps) if den > 0 else 1.0\n                d_ggo  = dice_np(pr[0].astype(bool), gt[0].astype(bool))\n                d_cons = dice_np(pr[1].astype(bool), gt[1].astype(bool))\n\n                # Overlays (verde = GGO, azul = Cons)\n                img8 = (img*255).clip(0,255).astype(np.uint8)\n                gt_ov = cv2.cvtColor(img8, cv2.COLOR_GRAY2BGR)\n                pr_ov = gt_ov.copy()\n                ggo_gt, cons_gt = gt[0].astype(bool),  gt[1].astype(bool)\n                ggo_pr, cons_pr = pr[0].astype(bool),  pr[1].astype(bool)\n\n                gt_ov[ggo_gt]  = (0.6*gt_ov[ggo_gt]  + 0.4*np.array([0,255,0])).astype(np.uint8)\n                gt_ov[cons_gt] = (0.6*gt_ov[cons_gt] + 0.4*np.array([0,0,255])).astype(np.uint8)\n                pr_ov[ggo_pr]  = (0.6*pr_ov[ggo_pr]  + 0.4*np.array([0,255,0])).astype(np.uint8)\n                pr_ov[cons_pr] = (0.6*pr_ov[cons_pr] + 0.4*np.array([0,0,255])).astype(np.uint8)\n\n                # Panel 1x3\n                fig, axs = plt.subplots(1,3, figsize=(12,4))\n                axs[0].imshow(img, cmap='gray'); axs[0].set_title('Input'); axs[0].axis('off')\n                axs[1].imshow(gt_ov[..., ::-1]); axs[1].set_title('Ground Truth'); axs[1].axis('off')\n                axs[2].imshow(pr_ov[..., ::-1])\n                axs[2].set_title(f'Prediction\\nDice GGO {d_ggo:.2f} | Cons {d_cons:.2f}')\n                axs[2].axis('off')\n                plt.tight_layout()\n\n                # Guardar + mostrar\n                fig.savefig(os.path.join(out_dir, f\"val_example_{saved:02d}.png\"), dpi=150)\n                if show:\n                    plt.show()\n                plt.close(fig)\n\n                saved += 1\n\nrender_examples(val_dl, n_samples=12, thr=0.45, show=True)\nprint(\"Saved in /kaggle/working/figs\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-18T12:45:36.907974Z","iopub.execute_input":"2025-09-18T12:45:36.908275Z","iopub.status.idle":"2025-09-18T12:45:36.927194Z","shell.execute_reply.started":"2025-09-18T12:45:36.908252Z","shell.execute_reply":"2025-09-18T12:45:36.926548Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#%%skip_if True\n\nimport pandas as pd\n\n@torch.no_grad()\ndef per_image_report(val_loader, thr=0.45, eps=1e-7, count_empty_as=1.0):\n    rows = []\n    model.eval()\n    for xb, yb in val_loader:\n        xb = xb.to(device); yb = yb.to(device)\n        logits = model(xb)\n        probs = torch.sigmoid(logits)\n        preds = (probs > thr).float()\n        B = xb.size(0)\n        for i in range(B):\n            row = {\"thr\": thr}\n            for c, cname in enumerate([\"ggo\", \"cons\"]):\n                p = preds[i,c]; t = yb[i,c]\n                inter = (p*t).sum().item()\n                den   = (p + t).sum().item()\n                if den == 0:\n                    dice = 1.0 if count_empty_as is not None else np.nan\n                    iou  = 1.0 if count_empty_as is not None else np.nan\n                else:\n                    dice = (2*inter + eps)/(den + eps)\n                    iou  = dice/(2 - dice)\n                row[f\"dice_{cname}\"] = dice\n                row[f\"iou_{cname}\"]  = iou\n            row[\"dice_mean\"] = np.nanmean([row[\"dice_ggo\"], row[\"dice_cons\"]])\n            row[\"iou_mean\"]  = np.nanmean([row[\"iou_ggo\"], row[\"iou_cons\"]])\n            rows.append(row)\n    return pd.DataFrame(rows)\n\ndf_report = per_image_report(val_dl, thr=0.45, count_empty_as=1.0)\ndf_report.to_csv(\"/kaggle/working/val_metrics_per_image.csv\", index=False)\n\nsummary = {\n    \"dice_ggo_mean\":  float(df_report[\"dice_ggo\"].mean()),\n    \"dice_cons_mean\": float(df_report[\"dice_cons\"].mean()),\n    \"dice_mean\":      float(df_report[\"dice_mean\"].mean()),\n    \"iou_ggo_mean\":   float(df_report[\"iou_ggo\"].mean()),\n    \"iou_cons_mean\":  float(df_report[\"iou_cons\"].mean()),\n    \"iou_mean\":       float(df_report[\"iou_mean\"].mean()),\n    \"n_images\":       int(len(df_report)),\n    \"thr\":            0.45\n}\nprint(\"Summary:\", summary)\nprint(\"CSV:\", \"/kaggle/working/val_metrics_per_image.csv\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-18T12:45:36.927786Z","iopub.execute_input":"2025-09-18T12:45:36.928028Z","iopub.status.idle":"2025-09-18T12:45:36.946631Z","shell.execute_reply.started":"2025-09-18T12:45:36.928005Z","shell.execute_reply":"2025-09-18T12:45:36.946051Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 6. Conclusions","metadata":{}},{"cell_type":"code","source":"# Mini viewer (2x2) para: Accuracy, Confusion, Loss, Results\nimport os\nimport matplotlib.pyplot as plt\nimport matplotlib.image as mpimg\n\nROOT = \"/kaggle/input/result-training-images\"\n\n# Orden fijo y candidatos por título\nentries = [\n    (\"Accuracy\",  [\"accuarcy.png\", \"accuracy.png\"]),  # intenta ambas\n    (\"Confusion\", [\"confusion.png\"]),\n    (\"Loss\",      [\"loss.png\"]),\n    (\"Results\",   [\"results.png\"]),\n]\n\ndef first_existing(root, candidates):\n    for name in candidates:\n        p = os.path.join(root, name)\n        if os.path.exists(p) and os.path.getsize(p) > 0:\n            return p\n    return None\n\n# Figura 2x2\nfig, axes = plt.subplots(2, 2, figsize=(12, 10))\naxes = axes.ravel()\n\nfor ax, (title, cands) in zip(axes, entries):\n    ax.axis(\"off\")\n    path = first_existing(ROOT, cands)\n    if path:\n        ax.imshow(mpimg.imread(path))\n        ax.set_title(title)\n    else:\n        ax.text(0.5, 0.5, f\"{title}\\n(file not found)\", ha=\"center\", va=\"center\", fontsize=12)\n\nfig.suptitle(\"Training Metrics Snapshots\", y=0.98)\nplt.tight_layout(rect=[0, 0, 1, 0.96])\n\nout_path = \"/kaggle/working/training_images_grid.png\"\nplt.savefig(out_path, dpi=150)\nprint(f\"[OK] Saved grid to {out_path}\")\n\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-19T04:24:28.021401Z","iopub.execute_input":"2025-09-19T04:24:28.021716Z","iopub.status.idle":"2025-09-19T04:24:29.630654Z","shell.execute_reply.started":"2025-09-19T04:24:28.021691Z","shell.execute_reply":"2025-09-19T04:24:29.629598Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 6.1 Accuarcy (per pixel)","metadata":{}},{"cell_type":"markdown","source":"The accuracy, although 99%, does not indicate good model results. Although the background and lung classes are not being targeted by the model, not classifying them already artificially increases their magnitude. If the model classified the entire image as background, it would have more than 95% accuracy.\nThe Grounded-glass opacity class reaches a coefficient (F1) of ~0.83, indicating that the vast majority of predicted pixels do match the actual mask, but its IoU of ~0.64 indicates that it still has many false values, making it difficult to delimit the edges.\nOn the other hand, the consolidation class is the minority in the image, so it does not have enough weight for correct identification, demonstrated by its coefficient of ~0.51 where only half of the predicted pixels are correct, and an IoU of ~0.36 confirms that the model is very conservative when it comes to masking a pixel with this class.\n\nAll of this, combined with the difference between the loss functions of the training and validation sets, indicates an overfitting problem with moderate bias and average variance.","metadata":{}},{"cell_type":"markdown","source":"### 6.2 Confusion Matrix","metadata":{}},{"cell_type":"markdown","source":"GGO\n\n- TN=38,712,350, FP=217,945, FN=386,092, TP=1,315,933\n\n- Precision 0.858, Recall 0.773, F1/Dice 0.813, Acc 0.985\n\nGood balance, some misses (FN) but many correct positives (TP).\n\nConsolidation\n\n- TN=40,377,709, FP=91,471, FN=74,901, TP=88,239\n\n- Precision 0.491, Recall 0.541, F1/Dice 0.515, Acc 0.996\n\nModel is conservative (low FP) but still misses many CONS (FN comparable to TP). High accuracy is again mostly background.","metadata":{}},{"cell_type":"markdown","source":"### 6.3 Loss train vs val","metadata":{}},{"cell_type":"markdown","source":"- Train loss drops smoothly as it progress 1.15 → 0.37.\n\n- Val loss improves early but then plateaus 0.90→0.76 with small bumps.\n\n- Persistent gap (≈ 0.39) ⇒ overfitting: the model probably memorizes train patterns better than it generalizes.\n\n- Notice that from the 30th-35th epoch the results improve vert smoothly or keeps between a range.","metadata":{}},{"cell_type":"markdown","source":"### 6.5 Underfitting/Overfitting?\n\nUnder/overfitting verdict\n\nThere is no evidence of underfitting respladed by our results, low train loss and Dice 0.70 however probably we have a small or mild overfitting, this due loss validation value do not improve as in training, also a stable Dice value but smaller than train. That behaviour, Validation accuracy near-perfect mainly due to huge TN counts (background), while foreground performance—especially CONS—lags (Dice ~0.52) are signals of overfitting.\n\n\n","metadata":{}},{"cell_type":"markdown","source":"## 7. Adjust hyperparameters and our model vs other\n\nTo amplify our vision we decide take action in two steps: \n\n1) Modify the hyperparameters of our current implementation.\n2) Use another model to compare the results.\n\nNote: This results are just to get more information about the performance of our current implementation and is not used it in test/submission.","metadata":{}},{"cell_type":"markdown","source":"## 7.1 Adjust hyperparameters\n\n### Hyperparameter changes — what and why\n\n- Loss variants:\n\n    - BCEWithLogits + Soft Dice → balances calibration (BCE) with overlap quality (Dice); stabilizes training and optimizes Dice directly.\n\n    - Weighted BCE+Dice (pos_weight) → counteracts class imbalance by up-weighting positive pixels per class.\n\n    - Tversky (α=0.7, β=0.3) → trades off FP vs FN; useful for small/rare lesions.\n\n- Upsampling head: Bilinear ↔ Transposed Conv → test whether learned upsampling sharpens boundaries vs smoother bilinear.\n\n- Learning rate: trial 1e-3 → 5e-4 on AdamW (wd=1e-4) → reduces oscillations/overfit; still decays with cosine annealing.\n\n- Decision threshold: tune around 0.45 → better Dice/accuracy trade-off for this dataset.\n\n- Epoch budget for ablations: ~20–30 epochs → quick feedback; keep full training for the main model.\n\nEverything else held constant: same data splits, augmentations, batch size, metrics—so gains are attributable to the tweak.","metadata":{}},{"cell_type":"code","source":"import time, copy, math\nimport numpy as np\nimport pandas as pd\nimport torch\nimport torch.nn as nn\n\n# ---------- Extra losses ----------\nclass BCEDiceLossWeighted(nn.Module):\n    \"\"\"BCEWithLogits (with pos_weight) + Soft Dice.\"\"\"\n    def __init__(self, pos_weight=None, eps=1e-7):\n        super().__init__()\n        self.bce = nn.BCEWithLogitsLoss(pos_weight=pos_weight)\n        self.eps = eps\n    def dice_soft(self, logits, targets):\n        probs = torch.sigmoid(logits)\n        inter = (probs * targets).sum(dim=(0,2,3))\n        den   = (probs + targets).sum(dim=(0,2,3)).clamp_min(self.eps)\n        dice  = (2*inter) / den\n        return 1 - dice.mean()\n    def forward(self, logits, targets):\n        return self.bce(logits, targets) + self.dice_soft(logits, targets)\n\nclass TverskyLoss(nn.Module):\n    \"\"\"Tversky loss (generalized Dice): alpha penalizes FN, beta penalizes FP.\"\"\"\n    def __init__(self, alpha=0.7, beta=0.3, eps=1e-7):\n        super().__init__()\n        self.alpha, self.beta, self.eps = float(alpha), float(beta), float(eps)\n    def forward(self, logits, targets):\n        probs = torch.sigmoid(logits)\n        TP = (probs * targets).sum(dim=(0,2,3))\n        FP = (probs * (1 - targets)).sum(dim=(0,2,3))\n        FN = ((1 - probs) * targets).sum(dim=(0,2,3))\n        tversky = (TP + self.eps) / (TP + self.alpha*FN + self.beta*FP + self.eps)\n        return 1 - tversky.mean()\n\n# ---------- Class imbalance → pos_weight for BCE ----------\ndef _compute_pos_weight_from_Y(Y_np):\n    # Y_np: (N,H,W,2) uint8\n    tot = Y_np.shape[0] * Y_np.shape[1] * Y_np.shape[2]\n    pos = Y_np.reshape(-1, 2).sum(axis=0).astype(np.float64)\n    pos = np.maximum(pos, 1.0)  # avoid div by zero\n    neg = tot - pos\n    pw = neg / pos  # per-channel pos_weight\n    return torch.tensor(pw, dtype=torch.float32, device=device)\n\n_pos_weight = _compute_pos_weight_from_Y(Y_train)\n\n# ---------- Local training helper (doesn't touch your globals) ----------\ndef train_one_model(build_model_fn, criterion, lr=1e-3, wd=1e-4, epochs=20,\n                    name=\"exp\", thr=None, use_cosine=True):\n    model_local = build_model_fn().to(device)\n    opt = torch.optim.AdamW(model_local.parameters(), lr=lr, weight_decay=wd)\n    sch = torch.optim.lr_scheduler.CosineAnnealingLR(opt, T_max=epochs) if use_cosine else None\n\n    best_val = -1.0\n    best_state = None\n    hist = {\"epoch\": [], \"train_dice\": [], \"val_dice\": [], \"train_loss\": [], \"val_loss\": []}\n    run_thr = float(thr if thr is not None else globals().get(\"THR\", 0.45))\n\n    # eval train set w/o augs (reuse if exists)\n    _train_eval_ds = NumpySegDataset(X_train, Y_train, train=False, apply_aug_fn=None)\n    _train_eval_dl = DataLoader(_train_eval_ds, batch_size=globals().get(\"BATCH\", 8),\n                                shuffle=False, num_workers=2, pin_memory=USE_CUDA)\n\n    start = time.time()\n    for ep in range(1, epochs+1):\n        model_local.train()\n        tr_loss_sum, n_items = 0.0, 0\n        for xb, yb in train_dl:\n            xb = xb.to(device, non_blocking=USE_CUDA); yb = yb.to(device, non_blocking=USE_CUDA)\n            opt.zero_grad(set_to_none=True)\n            with autocast_ctx():\n                logits = model_local(xb)\n                loss = criterion(logits, yb)\n            if USE_CUDA and scaler is not None:\n                scaler.scale(loss).backward()\n                scaler.step(opt); scaler.update()\n            else:\n                loss.backward(); opt.step()\n            tr_loss_sum += loss.item() * xb.size(0)\n            n_items += xb.size(0)\n\n        # eval\n        tr_loss_eval, tr_dice, _, _ = evaluate_epoch(model_local, _train_eval_dl, criterion, thr=run_thr)\n        va_loss_eval, va_dice, _, _ = evaluate_epoch(model_local, val_dl,        criterion, thr=run_thr)\n        if sch is not None: sch.step()\n\n        if va_dice > best_val:\n            best_val = va_dice\n            best_state = copy.deepcopy(model_local.state_dict())\n\n        hist[\"epoch\"].append(ep)\n        hist[\"train_loss\"].append(tr_loss_sum / max(n_items, 1))\n        hist[\"val_loss\"].append(va_loss_eval)\n        hist[\"train_dice\"].append(tr_dice)\n        hist[\"val_dice\"].append(va_dice)\n\n        print(f\"[{name}] Ep {ep:03d}/{epochs} | lr={opt.param_groups[0]['lr']:.2e} | \"\n              f\"tr_loss={hist['train_loss'][-1]:.4f} val_loss={va_loss_eval:.4f} \"\n              f\"tr_dice={tr_dice:.4f} val_dice={va_dice:.4f} best={best_val:.4f}\")\n\n    dur = time.time() - start\n    ckpt = f\"/kaggle/working/{name}.pt\"\n    if best_state is not None:\n        torch.save(best_state, ckpt)\n    return {\"name\": name, \"best_val_dice\": best_val, \"epochs\": epochs, \"lr\": lr,\n            \"wd\": wd, \"thr\": run_thr, \"seconds\": round(dur, 1), \"ckpt\": ckpt}, hist\n\ndef build_unet(bilinear=True):\n    return UNet2D(n_channels=1, n_classes=2, bilinear=bilinear)\n\nHP_EPOCHS = min(30, int(globals().get(\"EPOCHS\", 60)//2))  # short runs\nconfigs = [\n    {\"name\": \"unet_baseline_lr5e-4\",     \"bilinear\": True,  \"loss\": \"bce_dice\",          \"lr\": 5e-4, \"wd\": 1e-4},\n    {\"name\": \"unet_tversky_a0.7_b0.3\",   \"bilinear\": True,  \"loss\": \"tversky\",           \"lr\": 1e-3, \"wd\": 1e-4},\n    {\"name\": \"unet_weighted_bce_dice\",   \"bilinear\": True,  \"loss\": \"bce_dice_weighted\", \"lr\": 1e-3, \"wd\": 1e-4},\n    {\"name\": \"unet_deconv_upsample\",     \"bilinear\": False, \"loss\": \"bce_dice\",          \"lr\": 1e-3, \"wd\": 1e-4},\n]\n\nCONFIGS_TO_RUN = [c[\"name\"] for c in configs]  # or e.g. [\"unet_tversky_a0.7_b0.3\"]\n\nxp_rows, xp_hists = [], {}\nfor cfg in configs:\n    if cfg[\"name\"] not in CONFIGS_TO_RUN:\n        continue\n    # pick criterion\n    if cfg[\"loss\"] == \"tversky\":\n        criterion_hp = TverskyLoss(alpha=0.7, beta=0.3)\n    elif cfg[\"loss\"] == \"bce_dice_weighted\":\n        criterion_hp = BCEDiceLossWeighted(pos_weight=_pos_weight)\n    else:\n        criterion_hp = BCEDiceLoss()  # your original\n\n    res, hist = train_one_model(\n        build_model_fn=lambda: build_unet(bilinear=cfg[\"bilinear\"]),\n        criterion=criterion_hp, lr=cfg[\"lr\"], wd=cfg[\"wd\"],\n        epochs=HP_EPOCHS, name=cfg[\"name\"], thr=globals().get(\"THR\", 0.45)\n    )\n    xp_rows.append(res); xp_hists[cfg[\"name\"]] = hist\n\n# Save & display summary\nif xp_rows:\n    xp_df = pd.DataFrame(xp_rows).sort_values(\"best_val_dice\", ascending=False)\n    xp_csv = \"/kaggle/working/hparam_experiments.csv\"\n    xp_df.to_csv(xp_csv, index=False)\n    print(\"\\n=== Hyperparam Experiments (summary) ===\")\n    display(xp_df)\n    print(f\"Saved: {xp_csv}\")\nelse:\n    print(\"No experiments were run. Edit CONFIGS_TO_RUN.\")\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 7.2 DeepLab V3 vs our U-Net\n\n### But why DeepLab?\n\n- It brings multi-scale context via ASPP and dilated convolutions, which helps with lesions that vary from tiny GGOs to larger consolidations.\n\n- Provides a strong, widely-used segmentation baseline with a different inductive bias than U-Net (context-heavy vs skip-heavy), giving a fair test of whether global context is your current bottleneck.\n\n- Low engineering overhead (available in torchvision) and apples-to-apples evaluation (same data, loss, and metrics).\n\n- If DeepLab outperforms, it signals you might add context modules/receptive-field expansions to U-Net; if it ties, the task may be data-limited and further gains come from ETL/augmentation.\n\n### DeepLabV3 (comparison) — hyperparameter/config list\n\n- Backbone: ResNet-50 (torchvision DeepLabV3), aux_loss=False.\n\n- 0 Input: grayscale repeated to 3 channels.\n\n- Output: num_classes=2.\n\n- Loss: Weighted BCE+Dice (and Tversky as an alternative ablation).\n\n- Optimizer: AdamW, lr=3e-4, weight_decay=1e-4.\n\n- LR schedule: CosineAnnealingLR with T_max = epochs.\n\n- Epochs: up to 40 (or a short run for quick comparison).\n\n- Batch size: 8 (same as our U-Net).\n\n- Threshold for metrics: 0.45.\n\n- Augmentations & loaders: reuse exactly the same as U-Net.\n\n- Checkpoint: save best Val Dice to /kaggle/working/best_deeplabv3.pt.","metadata":{}},{"cell_type":"code","source":"# ============================================================\n# - Adapts 1-channel input by repeating to 3\n# - Sets num_classes=2\n# - Uses metrics/eval helpers\n# - Saves /kaggle/working/best_deeplabv3.pt\n# ============================================================\nimport torch\nimport torch.nn as nn\n\ntry:\n    from torchvision.models.segmentation import deeplabv3_resnet50\n    _has_tv = True\nexcept Exception as e:\n    print(\"torchvision not available for DeepLabV3:\", e)\n    _has_tv = False\n\nclass DeepLabV3Wrap(nn.Module):\n    \"\"\"Wraps torchvision DeepLabV3 to:\n       - Accept (B,1,H,W) by repeating channel to 3\n       - Return logits tensor (B,2,H,W) directly\n    \"\"\"\n    def __init__(self, num_classes=2):\n        super().__init__()\n        if not _has_tv:\n            raise RuntimeError(\"torchvision.deeplabv3_resnet50 not available.\")\n\n        # Try modern API\n        try:\n            self.net = deeplabv3_resnet50(weights=None, num_classes=num_classes, aux_loss=False)\n        except TypeError:\n            # Older API fallback\n            self.net = deeplabv3_resnet50(pretrained=False, progress=True, aux_loss=False)\n            # Replace classifier final conv out_channels to 2\n            if hasattr(self.net, \"classifier\") and isinstance(self.net.classifier, nn.Sequential):\n                # last module should be Conv2d(256, 21, 1)\n                for m in reversed(self.net.classifier):\n                    if isinstance(m, nn.Conv2d):\n                        if m.out_channels != num_classes:\n                            new_last = nn.Conv2d(m.in_channels, num_classes, kernel_size=1)\n                            m.in_channels, m.out_channels = new_last.in_channels, new_last.out_channels\n                            # replace in seq\n                            idx = list(self.net.classifier).index(m)\n                            self.net.classifier[idx] = new_last\n                        break\n\n    def forward(self, x):\n        # x: (B,1,H,W) → (B,3,H,W)\n        if x.size(1) == 1:\n            x = x.repeat(1, 3, 1, 1)\n        out = self.net(x)\n        # return raw logits tensor (torchvision returns dict in modern versions)\n        if isinstance(out, dict) and \"out\" in out:\n            return out[\"out\"]\n        return out  # assume tensor\n\ndef train_deeplab(epochs=None, lr=3e-4, wd=1e-4, thr=None, loss_name=\"bce_dice_weighted\"):\n    if not _has_tv:\n        print(\"Skipping DeepLabV3: torchvision unavailable.\")\n        return None, None\n\n    dlv3 = DeepLabV3Wrap(num_classes=2).to(device)\n\n    if loss_name == \"tversky\":\n        criterion_dl = TverskyLoss(alpha=0.7, beta=0.3)\n    elif loss_name == \"bce_dice_weighted\":\n        criterion_dl = BCEDiceLossWeighted(pos_weight=_pos_weight)\n    else:\n        criterion_dl = BCEDiceLoss()\n\n    opt = torch.optim.AdamW(dlv3.parameters(), lr=lr, weight_decay=wd)\n    ep = int(epochs if epochs is not None else min(40, int(globals().get(\"EPOCHS\", 60))))\n    run_thr = float(thr if thr is not None else globals().get(\"THR\", 0.45))\n    sch = torch.optim.lr_scheduler.CosineAnnealingLR(opt, T_max=ep)\n\n    best_val, best_state = -1.0, None\n    hist = {\"epoch\": [], \"train_dice\": [], \"val_dice\": [], \"train_loss\": [], \"val_loss\": []}\n\n    # train-eval loader for TRAIN w/o augs\n    _train_eval_ds = NumpySegDataset(X_train, Y_train, train=False, apply_aug_fn=None)\n    _train_eval_dl = DataLoader(_train_eval_ds, batch_size=globals().get(\"BATCH\", 8),\n                                shuffle=False, num_workers=2, pin_memory=USE_CUDA)\n\n    for e in range(1, ep+1):\n        dlv3.train()\n        tr_loss_sum, n_items = 0.0, 0\n        for xb, yb in train_dl:  # reuse your existing train_dl with augs\n            xb = xb.to(device, non_blocking=USE_CUDA); yb = yb.to(device, non_blocking=USE_CUDA)\n            opt.zero_grad(set_to_none=True)\n            with autocast_ctx():\n                logits = dlv3(xb)\n                loss = criterion_dl(logits, yb)\n            if USE_CUDA and scaler is not None:\n                scaler.scale(loss).backward()\n                scaler.step(opt); scaler.update()\n            else:\n                loss.backward(); opt.step()\n            tr_loss_sum += loss.item() * xb.size(0)\n            n_items += xb.size(0)\n\n        tr_loss_eval, tr_dice, _, _ = evaluate_epoch(dlv3, _train_eval_dl, criterion_dl, thr=run_thr)\n        va_loss_eval, va_dice, _, _ = evaluate_epoch(dlv3, val_dl,        criterion_dl, thr=run_thr)\n        sch.step()\n\n        if va_dice > best_val:\n            best_val = va_dice\n            best_state = copy.deepcopy(dlv3.state_dict())\n\n        hist[\"epoch\"].append(e)\n        hist[\"train_loss\"].append(tr_loss_sum / max(n_items, 1))\n        hist[\"val_loss\"].append(va_loss_eval)\n        hist[\"train_dice\"].append(tr_dice)\n        hist[\"val_dice\"].append(va_dice)\n\n        print(f\"[deeplabv3] Ep {e:03d}/{ep} | lr={opt.param_groups[0]['lr']:.2e} | \"\n              f\"tr_loss={hist['train_loss'][-1]:.4f} val_loss={va_loss_eval:.4f} \"\n              f\"tr_dice={tr_dice:.4f} val_dice={va_dice:.4f} best={best_val:.4f}\")\n\n    ckpt = \"/kaggle/working/best_deeplabv3.pt\"\n    if best_state is not None:\n        torch.save(best_state, ckpt)\n    return {\"name\": \"deeplabv3_resnet50\", \"best_val_dice\": best_val, \"epochs\": ep,\n            \"lr\": lr, \"wd\": wd, \"thr\": run_thr, \"ckpt\": ckpt}, hist\n\n# ----  DeepLabV3 comparison ----\ndl_summary, dl_hist = train_deeplab(epochs=min(40, int(globals().get(\"EPOCHS\", 60))),\n                                    lr=3e-4, wd=1e-4, loss_name=\"bce_dice_weighted\")\n\n# ---- Compare U-Net best vs DeepLabV3 ----\nrows = []\ntry:\n    rows.append({\"model\": \"unet_yours_best\", \"best_val_dice\": float(globals().get(\"best_val\", -1.0)),\n                 \"ckpt\": \"/kaggle/working/best_unet.pt\"})\nexcept Exception:\n    pass\nif dl_summary is not None:\n    rows.append({\"model\": dl_summary[\"name\"], \"best_val_dice\": dl_summary[\"best_val_dice\"],\n                 \"ckpt\": dl_summary[\"ckpt\"]})\nif rows:\n    cmp_df = pd.DataFrame(rows).sort_values(\"best_val_dice\", ascending=False)\n    print(\"\\n=== Model comparison (val Dice) ===\")\n    display(cmp_df)\n    cmp_df.to_csv(\"/kaggle/working/model_comparison.csv\", index=False)\n    print(\"Saved: /kaggle/working/model_comparison.csv\")\nelse:\n    print(\"No models to compare.\")\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 8. Overall Conclussions\n\n\nOverall and according with our results we get a more than acceptable results, and we detect the following points to support our good performance at our first baseline implementeation:\n\n- Stable training: no divergence; train loss ↓ 1.15 → 0.37, val metrics trend upward and stabilize.\n\n- Meaningful segmentation: not an all-background cheat—your val Dicē ≈ 0.664 with GGO Dice ≈ 0.813 shows the model really finds lesions.\n\n- Gaps are modest: train vs val Dice gap ~4 pp = mild overfitting, not catastrophic.\n\n- Confusions make sense: high TNs (class imbalance), reasonable GGO TP/FP tradeoff; CONS is harder (typical).\n\nFor our purpose we unlock a significant archivement ss a first baseline, results are more than acceptable and provide a solid platform for targeted improvements—especially for CONS and generalization.\n\n### Pros \n\nOur model has strengths that will let us put hands on imprvement.\n\n- Comes from U-Net base implementation so is not a very technical or difficult variation.\n\n- Convergent and stable training with a sensible schedule (AdamW + cosine).\n\n- Strong GGO segmentation (Dice ~0.81) with balanced precision/recall.\n\n- Reproducibility: fixed seed, deterministic cuDNN, clear artifact logging (CSV + figures + best checkpoint).\n\n- Operating point already serviceable: a single global threshold (τ=0.45) yields decent mean Dice without heavy tuning.\n\n- Good specificity / low FP consistent with medical expectations.\n\n### Cons\nBut also there is prove that mild overfitting and in consecuence, oportunity areas of opportunity like:\n\n- Mild overfitting: persistent train–val gaps; val loss plateaus around 0.76.\n\n- CONS underperformance: small and low-contrast regions yield lower Dice/recall than GGO.\n\n- Accuracy is misleading: dominated by background TNs; Dice/IoU and balanced accuracy are more informative.\n\n- Augmentation & sampling not class-aware: generic aug helps, but without lesion-aware sampling, CONS remains rare in batches.\n\n### Future work and improvements\n\n- Openess to feedback: First of all we are opened to receibe feedback from experts or more experienced people in the filed as our professors, community, etc... this in order to expand our vision and improve our model.\n\n- Loss tilted toward FN: Try Tversky/Focal-Tversky with α=0.3, β=0.7 (β>α penalizes FN). Or keep BCE+Dice but set pos_weight>1 only for CONS.\n\n- Early stopping on val Dice.\n\n- Try Tversky/Focal-Tversky with α=0.3, β=0.7 (β>α penalizes FN). Or keep BCE+Dice but set pos_weight>1 only for CONS.\n\n- Lesion-aware sampling: Ensure a fraction of training crops must contain CONS pixels (min-area filter or oversample CONS-positive slices).\n\n- Better report diagnosis: Alongside pixel accuracy, include balanced accuracy, Cohen’s κ, per-image median Dice, and PR curves.\n\n- More ","metadata":{}},{"cell_type":"markdown","source":"## Test\n\nThis section is dedicated to run our model (U-Net) in Test data set in order to generate our results and submit it to the competiotion. ","metadata":{}},{"cell_type":"code","source":"# ============================\n# TEST → SUBMISSION (Id,Predicted)\n# ============================\nimport os, numpy as np, pandas as pd, torch\n\ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n\n# 1) Load weights if is necessary\nckpt_path = \"/kaggle/working/best_unet.pt\"\nif os.path.exists(ckpt_path):\n    state = torch.load(ckpt_path, map_location=device)\n    model.load_state_dict(state)\nmodel.eval()\n\ntry:\n    X_test  # (N,512,512,1) float32 [0,1]\nexcept NameError:\n    prep = \"/kaggle/working/prepared\"\n    X_test = np.load(os.path.join(prep, \"test_medseg.npz\"))[\"X\"]\n\n# 3) Prediction\nTHR   = float(globals().get(\"THR\", 0.45))\nBATCH = int(globals().get(\"BATCH\", 8))\n\n@torch.no_grad()\ndef predict_binary(model, X, thr=THR, batch=BATCH, device=device):\n    outs = []\n    for i in range(0, len(X), batch):\n        xb = torch.from_numpy(X[i:i+batch]).permute(0,3,1,2).float().to(device)  # (B,1,H,W)\n        logits = model(xb)                               # (B,2,H,W)\n        probs  = torch.sigmoid(logits).cpu().numpy()     # (B,2,H,W)\n        preds  = (probs >= thr).astype(np.uint8)         # binario\n        outs.append(preds)\n    return np.concatenate(outs, axis=0)                  # (N,2,H,W)\n\npreds = predict_binary(model, X_test)\nN, C, H, W = preds.shape\nassert C == 2, f\"Esperaba 2 clases, obtuve {C}\"\n\n# 4) Keeps the requested format (N,H,W,2) before ravel()\nif preds.shape == (N, C, H, W):                         # NCHW → NHWC\n    preds_nhwc = np.transpose(preds, (0, 2, 3, 1)).copy()\nelse:\n    preds_nhwc = preds\n\n# Sanity check \nu = np.unique(preds_nhwc)\nassert set(u.tolist()).issubset({0,1}), f\"No binary predictions: {u}\"\n\n# 5) Build CSV\nflat = preds_nhwc.ravel(order='C').astype(int)          # (N*H*W*2,)\nids  = np.arange(flat.size, dtype=int)\n\nsub_df = pd.DataFrame({\"Id\": ids, \"Predicted\": flat}).set_index(\"Id\")\nout_path = \"/kaggle/working/submission.csv\"\nsub_df.to_csv(out_path)\n\nprint(f\"[OK] submission.csv guardado en {out_path}\")\nprint(\"Expected size:\", N*H*W*2, \"| rows:\", len(sub_df))\nprint(sub_df.head())\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## References\n\nO. Ronneberger, P. Fischer, and T. Brox, “U-Net: Convolutional Networks for Biomedical Image Segmentation,” in Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015, LNCS, vol. 9351, N. Navab, J. Hornegger, W. M. Wells, and A. F. Frangi, Eds. Cham: Springer, 2015, pp. 234–241. doi: 10.1007/978-3-319-24574-4_28.\n\nChen, L.-C., Papandreou, G., Kokkinos, I., Murphy, K., & Yuille, A. L. (2017). DeepLab: Semantic image segmentation with deep convolutional nets, atrous convolution, and fully connected CRFs (Version 2). arXiv. https://doi.org/10.48550/arXiv.1606.00915\n\nMaftouni M. (2020). PyTorch Baseline for Semantic Segmentation. Kaggle. https://www.kaggle.com/code/maedemaftouni/pytorch-baseline-for-semantic-segmentation","metadata":{}}]}