{"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":[{"sourceType":"competition","sourceId":99552,"databundleVersionId":13441085}],"dockerImageVersionId":31089,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"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","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-09-02T14:18:35.169483Z","iopub.execute_input":"2025-09-02T14:18:35.169807Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ============================================================\n# RSNA Intracranial Aneurysm Detection: Geometry metrics from\n# COW NIfTI segmentations + stats by named artery\n# ============================================================\n\n# ---------- Cell 0: Install lightweight deps ----------\n# (Kaggle usually has numpy/pandas/scipy; nibabel/tqdm added here)\n!pip -q install nibabel tqdm\n\n# ---------- Cell 1: Imports & basic setup ----------\nimport os, glob\nfrom pathlib import Path\nimport numpy as np\nimport pandas as pd\nimport nibabel as nib\nfrom tqdm import tqdm\nfrom IPython.display import display, clear_output\n\n# Competition root and standard files\nCOMP_DIR   = \"/kaggle/input/rsna-intracranial-aneurysm-detection\"\nTRAIN_CSV  = os.path.join(COMP_DIR, \"train.csv\")   # Adjust if your train file lives elsewhere\n\n# Where to write outputs (persist with notebook versioning)\nWORK_DIR   = \"/kaggle/working\"\nOUT_DIR    = os.path.join(WORK_DIR, \"aneurysm_metrics_fast\")\nOUT_CSV    = os.path.join(OUT_DIR, \"segment_metrics.csv\")\n\n# ------------- Segmentations -------------\n# We auto-discover any *_cowseg.nii* anywhere under /kaggle/input.\n# If you prefer a fixed folder, set SEG_ROOT and use FILE_GLOB under it.\nSEG_ROOT   = \"/kaggle/input\"           # search base\nFILE_GLOB  = \"**/*_cowseg.nii*\"        # pattern for labeled segmentations\n\n# ------------- Metrics controls -------------\nMIN_VOL_MM3 = 5.0     # ignore tiny blobs\nSHOW_EVERY  = 1       # how often to refresh console & preview df\n\nos.makedirs(OUT_DIR, exist_ok=True)\n\n# ---------- Cell 2: NIfTI I/O and fast surface area ----------\ndef nii_read(path):\n    \"\"\"\n    Load a NIfTI and return (data_uint8, spacing_zyx_mm)\n    \"\"\"\n    nii = nib.load(path)\n    data = np.asarray(nii.get_fdata(), dtype=np.uint8)\n    spacing = tuple(float(z) for z in nii.header.get_zooms()[:3])  # (z, y, x) in mm\n    return data, spacing\n\ndef fast_surface_area(binary, spacing):\n    \"\"\"\n    6-neighborhood exposed-face count * face areas.\n    Approximates surface area quickly without marching cubes.\n    \"\"\"\n    if binary.dtype != bool:\n        binary = binary.astype(bool)\n    dz, dy, dx = spacing\n    ax = dy * dz  # faces normal to X\n    ay = dx * dz  # faces normal to Y\n    az = dx * dy  # faces normal to Z\n\n    exposed = 0.0\n    # +X\n    nb = np.zeros_like(binary); nb[..., :-1] = binary[..., 1:]\n    exposed += ax * np.sum(binary & ~nb)\n    # -X\n    nb = np.zeros_like(binary); nb[..., 1:] = binary[..., :-1]\n    exposed += ax * np.sum(binary & ~nb)\n    # +Y\n    nb = np.zeros_like(binary); nb[:, :-1, :] = binary[:, 1:, :]\n    exposed += ay * np.sum(binary & ~nb)\n    # -Y\n    nb = np.zeros_like(binary); nb[:, 1:, :] = binary[:, :-1, :]\n    exposed += ay * np.sum(binary & ~nb)\n    # +Z\n    nb = np.zeros_like(binary); nb[:-1, :, :] = binary[1:, :, :]\n    exposed += az * np.sum(binary & ~nb)\n    # -Z\n    nb = np.zeros_like(binary); nb[1:, :, :] = binary[:-1, :, :]\n    exposed += az * np.sum(binary & ~nb)\n\n    return float(exposed)\n\ndef compute_metrics_for_label(mask_bool, spacing):\n    \"\"\"\n    Basic 3D metrics for a labeled component:\n    - Volume (mm^3)\n    - Surface area approx. (mm^2) via exposed faces\n    - Max diameter approx. (mm): 2*max distance from centroid\n    - Elongation (PCA): sqrt(lambda1/lambda3)\n    - Sphericity: (pi^(1/3) * (6V)^(2/3)) / A\n    \"\"\"\n    voxvol = spacing[0] * spacing[1] * spacing[2]\n    nvox = int(mask_bool.sum())\n    volume_mm3 = nvox * voxvol\n    if nvox == 0:\n        return dict(volume_mm3=0.0, surface_mm2=0.0, max_diameter_mm=0.0,\n                    elongation=np.nan, sphericity=np.nan)\n\n    area_mm2 = fast_surface_area(mask_bool, spacing)\n\n    zz, yy, xx = np.where(mask_bool)\n    coords_mm = np.vstack([zz*spacing[0], yy*spacing[1], xx*spacing[2]]).T\n\n    center = coords_mm.mean(axis=0, keepdims=True)\n    dists = np.linalg.norm(coords_mm - center, axis=1)\n    max_diameter_mm = float(dists.max() * 2.0) if dists.size else 0.0\n\n    # PCA-based elongation\n    if coords_mm.shape[0] >= 5:\n        C = coords_mm - center\n        cov = (C.T @ C) / max(len(C)-1, 1)\n        w = np.linalg.eigvalsh(cov)\n        w = np.sort(np.clip(w, 1e-12, None))[::-1]  # lambda1>=lambda2>=lambda3>0\n        elongation = float(np.sqrt(w[0]/w[-1]))\n    else:\n        elongation = np.nan\n\n    # Sphericity\n    if area_mm2 > 0 and volume_mm3 > 0:\n        sphericity = float((np.pi**(1/3.0)) * ((6.0*volume_mm3)**(2/3.0)) / area_mm2)\n    else:\n        sphericity = np.nan\n\n    return dict(volume_mm3=float(volume_mm3),\n                surface_mm2=float(area_mm2),\n                max_diameter_mm=float(max_diameter_mm),\n                elongation=elongation,\n                sphericity=sphericity)\n\n# ---------- Cell 3: Discover NIfTI segmentation files ----------\npaths = sorted(glob.glob(os.path.join(SEG_ROOT, FILE_GLOB), recursive=True))\nprint(f\"Found {len(paths)} volumes matching pattern: {FILE_GLOB}\")\n\nif len(paths) == 0:\n    raise FileNotFoundError(\n        \"No NIfTI segmentations were found. \"\n        \"Make sure you've added a dataset with files named like *_cowseg.nii.gz under /kaggle/input/.\\n\"\n        \"You can also change SEG_ROOT/FILE_GLOB near the top of this notebook.\"\n    )\n\n# ---------- Cell 4: Iterate segmentations → per-label metrics CSV ----------\nif os.path.exists(OUT_CSV):\n    df_all = pd.read_csv(OUT_CSV)\n    processed = set((df_all[\"SeriesInstanceUID\"].astype(str) + \"_\" + df_all[\"Label\"].astype(str)).tolist())\nelse:\n    df_all = pd.DataFrame(columns=[\n        \"SeriesInstanceUID\",\"Label\",\n        \"Volume_mm3\",\"SurfaceArea_mm2\",\"MaxDiameter_mm\",\n        \"Elongation\",\"Sphericity\"\n    ])\n    processed = set()\n\nfor i, nii_path in enumerate(tqdm(paths, desc=\"Volumes\", unit=\"vol\")):\n    seg, spacing = nii_read(nii_path)\n    uid = Path(nii_path).stem.replace(\"_cowseg\",\"\")\n    if uid == Path(nii_path).stem:\n        # if suffix not present, fallback to folder name\n        uid = Path(nii_path).parent.name\n\n    labels = np.unique(seg)\n    labels = labels[labels != 0]\n    if labels.size == 0:\n        continue\n\n    new_rows = []\n    for lbl in labels:\n        key = f\"{uid}_{int(lbl)}\"\n        if key in processed:\n            continue\n        mask = (seg == lbl)\n\n        # quick volume filter\n        n_vox = int(mask.sum())\n        vol_mm3 = n_vox * spacing[0] * spacing[1] * spacing[2]\n        if vol_mm3 < MIN_VOL_MM3:\n            continue\n\n        metrics = compute_metrics_for_label(mask, spacing)\n        row = {\n            \"SeriesInstanceUID\": uid,\n            \"Label\": int(lbl),\n            \"Volume_mm3\": metrics[\"volume_mm3\"],\n            \"SurfaceArea_mm2\": metrics[\"surface_mm2\"],\n            \"MaxDiameter_mm\": metrics[\"max_diameter_mm\"],\n            \"Elongation\": metrics[\"elongation\"],\n            \"Sphericity\": metrics[\"sphericity\"],\n        }\n        new_rows.append(row)\n        processed.add(key)\n\n    if new_rows:\n        df_all = pd.concat([df_all, pd.DataFrame(new_rows)], ignore_index=True)\n        Path(OUT_CSV).parent.mkdir(parents=True, exist_ok=True)\n        df_all.to_csv(OUT_CSV, index=False)\n\n    if ((i+1) % max(1, SHOW_EVERY) == 0) or (i+1 == len(paths)):\n        clear_output(wait=True)\n        print(f\"Processed {i+1}/{len(paths)} volumes. Last UID: {uid}\")\n        if not df_all.empty:\n            display(df_all.tail(8))\n\nprint(\"Saved metrics CSV to:\", OUT_CSV)\n\n# ============================================================\n# Stats by artery name (merge with train.csv)\n# ============================================================\n\n# ---------- Cell 5: Imports for stats ----------\nfrom scipy import stats\n\n# Output folder for stats\nSTATS_OUT_DIR = os.path.join(WORK_DIR, \"segment_stats_out_by_name\")\nos.makedirs(STATS_OUT_DIR, exist_ok=True)\n\n# Label map (edit if your label indices → artery names differ)\nLABELS = {\n    1: \"Other Posterior Circulation\",\n    2: \"Basilar Tip\",\n    3: \"Right Posterior Communicating Artery\",\n    4: \"Left Posterior Communicating Artery\",\n    5: \"Right Infraclinoid Internal Carotid Artery\",\n    6: \"Left Infraclinoid Internal Carotid Artery\",\n    7: \"Right Supraclinoid Internal Carotid Artery\",\n    8: \"Left Supraclinoid Internal Carotid Artery\",\n    9: \"Right Middle Cerebral Artery\",\n    10: \"Left Middle Cerebral Artery\",\n    11: \"Right Anterior Cerebral Artery\",\n    12: \"Left Anterior Cerebral Artery\",\n    13: \"Anterior Communicating Artery\",\n}\nFEATURES = ['Volume_mm3', 'SurfaceArea_mm2', 'MaxDiameter_mm', 'Elongation', 'Sphericity']\n\nMIN_SAMPLES_PER_GROUP = 5  # per group, per-label\nALPHA = 0.05\nEFFECT_MIN = 0.147  # small but noticeable rank-biserial correlation\n\n# ---------- Cell 6: Load CSVs and harmonize UID column ----------\nmetrics = pd.read_csv(OUT_CSV)\n\nif not os.path.exists(TRAIN_CSV):\n    raise FileNotFoundError(\n        f\"Could not find train.csv at: {TRAIN_CSV}\\n\"\n        \"If your train CSV is elsewhere, update TRAIN_CSV near the top.\"\n    )\n\ntrain = pd.read_csv(TRAIN_CSV)\n\n# Try to find a column equivalent to SeriesInstanceUID\nuid_candidates = [\n    \"SeriesInstanceUID\", \"StudyInstanceUID\", \"StudyUID\",\n    \"series_id\", \"study_id\", \"SeriesID\", \"StudyID\"\n]\nuid_found = None\nfor c in uid_candidates:\n    if c in train.columns:\n        uid_found = c\n        break\n\nif uid_found is None:\n    raise ValueError(\n        \"Could not find a UID column in train.csv. \"\n        f\"Tried: {uid_candidates}. Please rename your UID column to 'SeriesInstanceUID' \"\n        \"or extend uid_candidates above.\"\n    )\n\nif uid_found != \"SeriesInstanceUID\":\n    train = train.rename(columns={uid_found: \"SeriesInstanceUID\"})\n\n# Identify artery columns by exact name match to LABELS values\nartery_cols = [name for name in LABELS.values() if name in train.columns]\nif not artery_cols:\n    raise ValueError(\n        \"No artery-name columns were found in train.csv.\\n\"\n        \"Expected headers like those in LABELS (e.g., 'Left Middle Cerebral Artery').\\n\"\n        \"Please align LABELS values to your train.csv column names.\"\n    )\n\n# Optionally include a global 'Aneurysm Present' if available\ncols_to_merge = ['SeriesInstanceUID'] + artery_cols + ([c for c in ['Aneurysm Present'] if c in train.columns])\n\nmerged = metrics.merge(train[cols_to_merge], on='SeriesInstanceUID', how='left')\nmerged['LabelName'] = merged['Label'].map(LABELS)\n\n# Per-row flag: aneurysm present in THIS labeled artery\ndef row_flag(row, frame_cols):\n    ln = row['LabelName']\n    if ln in frame_cols:\n        val = row[ln]\n        try:\n            return int(val)\n        except Exception:\n            return 0\n    return np.nan  # column missing — drop later\n\nmerged['Aneurysm_In_This_Label'] = merged.apply(lambda r: row_flag(r, merged.columns), axis=1).astype('float')\n\n# Filter valid rows\nvalid = merged[merged['Aneurysm_In_This_Label'].notna()].copy()\nvalid['Aneurysm_In_This_Label'] = valid['Aneurysm_In_This_Label'].astype(int)\n\n# ---------- Cell 7: Mann–Whitney utilities ----------\ndef rank_biserial_from_u(u_stat, n1, n0):\n    return (u_stat / (n1 * n0) - 0.5) * 2\n\ndef mannwhitney_summary(group1, group0):\n    g1 = np.asarray(pd.Series(group1).dropna().values, dtype=float)\n    g0 = np.asarray(pd.Series(group0).dropna().values, dtype=float)\n    if len(g1)==0 or len(g0)==0:\n        return dict(n1=len(g1), n0=len(g0), median1=np.nan, median0=np.nan, U=np.nan, p=np.nan, rbc=np.nan)\n    U, p = stats.mannwhitneyu(g1, g0, alternative='two-sided')\n    rbc = rank_biserial_from_u(U, len(g1), len(g0))\n    return dict(n1=len(g1), n0=len(g0), median1=np.median(g1), median0=np.median(g0), U=U, p=p, rbc=rbc)\n\n# ---------- Cell 8: GLOBAL comparison (flag by specific LabelName) ----------\nglobal_rows = []\nfor feat in FEATURES:\n    res = mannwhitney_summary(\n        valid.loc[valid['Aneurysm_In_This_Label']==1, feat],\n        valid.loc[valid['Aneurysm_In_This_Label']==0, feat]\n    )\n    global_rows.append({\n        'Feature': feat,\n        'Median_Aneurysm(ThisLabel)': res['median1'],\n        'Median_Normal(ThisLabel)': res['median0'],\n        'U_stat': res['U'],\n        'p_value': res['p'],\n        'Rank_Biserial_Corr': res['rbc'],\n        'N_Aneurysm(ThisLabel)': res['n1'],\n        'N_Normal(ThisLabel)': res['n0'],\n    })\nglobal_df = pd.DataFrame(global_rows).sort_values('p_value')\nglobal_df.to_csv(Path(STATS_OUT_DIR)/\"global_by_labelname_stats.csv\", index=False)\n\n# ---------- Cell 9: Per-label (per artery id) comparison ----------\nperlabel_rows = []\nfor lid, lname in LABELS.items():\n    sub = valid[valid['Label'] == lid]\n    if sub.empty:\n        continue\n    g1 = sub[sub['Aneurysm_In_This_Label']==1]\n    g0 = sub[sub['Aneurysm_In_This_Label']==0]\n    if len(g1) < MIN_SAMPLES_PER_GROUP or len(g0) < MIN_SAMPLES_PER_GROUP:\n        continue\n    for feat in FEATURES:\n        res = mannwhitney_summary(g1[feat], g0[feat])\n        perlabel_rows.append({\n            'Label': lid,\n            'LabelName': lname,\n            'Feature': feat,\n            'Median_Aneurysm(ThisLabel)': res['median1'],\n            'Median_Normal(ThisLabel)': res['median0'],\n            'U_stat': res['U'],\n            'p_value': res['p'],\n            'Rank_Biserial_Corr': res['rbc'],\n            'N_Aneurysm(ThisLabel)': res['n1'],\n            'N_Normal(ThisLabel)': res['n0'],\n        })\nperlabel_df = pd.DataFrame(perlabel_rows).sort_values(['Label','p_value'])\nperlabel_df.to_csv(Path(STATS_OUT_DIR)/\"perlabel_by_labelname_stats.csv\", index=False)\n\n# ---------- Cell 10: Significance filters & console summary ----------\ndef mark_sig(df):\n    return (df['p_value'] < ALPHA) & (df['Rank_Biserial_Corr'].abs() >= EFFECT_MIN)\n\nsig_global = global_df[mark_sig(global_df)].copy()\nsig_perlabel = perlabel_df[mark_sig(perlabel_df)].copy()\n\nsig_global.to_csv(Path(STATS_OUT_DIR)/\"global_by_labelname_stats_significant.csv\", index=False)\nsig_perlabel.to_csv(Path(STATS_OUT_DIR)/\"perlabel_by_labelname_stats_significant.csv\", index=False)\n\nprint(\">> GLOBAL by-artery (flag specific to LabelName):\")\nprint(global_df.to_string(index=False))\nprint(\"\\nSignificant (p < %.3f & |r| >= %.3f):\" % (ALPHA, EFFECT_MIN))\nprint(sig_global.to_string(index=False) if not sig_global.empty else \"(none)\")\n\nprint(\"\\n>> PER-LABEL (only labels with n≥%d per group):\" % MIN_SAMPLES_PER_GROUP)\nif perlabel_df.empty:\n    print(\"(Not enough samples per label)\")\nelse:\n    print(perlabel_df.sort_values('p_value').head(12).to_string(index=False))\n    print(\"\\nSignificant (p < %.3f & |r| >= %.3f):\" % (ALPHA, EFFECT_MIN))\n    print(sig_perlabel.sort_values('p_value').to_string(index=False) if not sig_perlabel.empty else \"(none)\")\n\nprint(\"\\nMetrics CSV:\", Path(OUT_CSV).resolve())\nprint(\"Stats out dir:\", Path(STATS_OUT_DIR).resolve())\n","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}