{"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":"nvidiaTeslaT4","dataSources":[{"sourceId":36363,"databundleVersionId":4050810,"sourceType":"competition"},{"sourceId":71549,"databundleVersionId":8561470,"sourceType":"competition"}],"dockerImageVersionId":31090,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# === Cell 0: Imports, config, seeding, util print ===\nimport os, sys, glob, random, textwrap, traceback\nfrom pathlib import Path\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport matplotlib as mpl\nimport SimpleITK as sitk\nfrom scipy import ndimage as ndi\n\n# Matplotlib aesthetics\nmpl.rcParams[\"figure.dpi\"] = 140\nmpl.rcParams[\"font.size\"] = 11\n\ndef seed_all(s=1337):\n random.seed(s); np.random.seed(s)\nseed_all(1337)\n\nprint(\"Python:\", sys.version)\nprint(\"SimpleITK:\", sitk.Version_VersionString())\n\n# CT windowing helpers\ndef window_image(img_np, center=400, width=1800):\n lo = center - width/2.0\n hi = center + width/2.0\n x = np.clip(img_np, lo, hi)\n x = (x - lo) / (hi - lo + 1e-6)\n return x.astype(np.float32)\n\n# quick border extractor for label images\ndef boundary_from_label(lbl2d):\n # binary edge for each class != 0 (union), using morphological erosion\n mask = (lbl2d > 0).astype(np.uint8)\n er = ndi.binary_erosion(mask, structure=np.ones((3,3)), iterations=1)\n edge = (mask ^ er).astype(np.uint8)\n return edge\n\n# multi-color overlay of label map\ndef overlay_segmentation(img2d_0to1, lbl2d, alpha=0.35):\n h,w = img2d_0to1.shape\n base = np.stack([img2d_0to1]*3, axis=-1) # gray->RGB\n uniq = sorted([u for u in np.unique(lbl2d) if u>0])\n # stable color palette across runs\n rng = np.random.default_rng(42)\n colors = rng.random((max(uniq+[1])+1, 3)) # indexed by label id\n # emphasize C1..C7 labels (1..7) with a fixed palette\n fixed = {\n 1:(0.85,0.10,0.10), 2:(0.10,0.55,0.85), 3:(0.10,0.75,0.30),\n 4:(0.75,0.55,0.10), 5:(0.60,0.10,0.75), 6:(0.15,0.80,0.80), 7:(0.90,0.30,0.30)\n }\n for k,v in fixed.items(): colors[k]=v\n out = base.copy()\n for u in uniq:\n mask = (lbl2d==u)\n color = colors[u]\n out[mask] = (1-alpha)*out[mask] + alpha*np.array(color)\n # thin boundary lines for crispness\n edge = boundary_from_label(lbl2d)\n out[edge.astype(bool)] = (0,0,0) # black edge\n return out\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T21:44:51.389221Z","iopub.execute_input":"2025-09-24T21:44:51.38953Z","iopub.status.idle":"2025-09-24T21:44:54.963096Z","shell.execute_reply.started":"2025-09-24T21:44:51.389499Z","shell.execute_reply":"2025-09-24T21:44:54.96211Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Cell 0P: Patch overlay_segmentation to handle arbitrary label IDs safely ===\ndef overlay_segmentation(img2d_0to1, lbl2d, alpha=0.35):\n # img2d_0to1: (H,W) in [0,1], lbl2d: (H,W) integer labels\n h, w = img2d_0to1.shape\n base = np.stack([img2d_0to1]*3, axis=-1)\n lbl2d = lbl2d.astype(np.int32)\n\n uniq = [int(u) for u in np.unique(lbl2d) if u > 0]\n if not uniq:\n return base # no labels to overlay\n\n # Ensure palette big enough for cervical (1..7) and thoracic (8..19)\n max_label = max(max(uniq), 19) # allocate up to 19\n rng = np.random.default_rng(42)\n colors = rng.random((max_label + 1, 3))\n\n # Fixed readable colors for C1..C7\n fixed = {\n 1:(0.85,0.10,0.10), 2:(0.10,0.55,0.85), 3:(0.10,0.75,0.30),\n 4:(0.75,0.55,0.10), 5:(0.60,0.10,0.75), 6:(0.15,0.80,0.80),\n 7:(0.90,0.30,0.30)\n }\n for k, v in fixed.items():\n if k <= max_label:\n colors[k] = v\n\n out = base.copy()\n for u in uniq:\n mask = (lbl2d == u)\n out[mask] = (1 - alpha) * out[mask] + alpha * colors[u]\n\n # crisp boundary\n edge = boundary_from_label(lbl2d)\n out[edge.astype(bool)] = (0, 0, 0)\n return out\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T21:53:36.022217Z","iopub.execute_input":"2025-09-24T21:53:36.02309Z","iopub.status.idle":"2025-09-24T21:53:36.030933Z","shell.execute_reply.started":"2025-09-24T21:53:36.02306Z","shell.execute_reply":"2025-09-24T21:53:36.030032Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Cell 1: Discover RSNA 2022 CT dataset with vertebra segmentations ===\n\nCANDIDATE_ROOTS = [\n \"/kaggle/input/rsna-2022-cervical-spine-fracture-detection\",\n \"/kaggle/input/rsna-2022-cervical-spine-fracture-detection-2\", # safety\n \"/kaggle/input/rsna-2022-spine-fracture-detection\", # occasional alias\n]\n\nDATA_ROOT = None\nfor r in CANDIDATE_ROOTS:\n if Path(r).exists():\n DATA_ROOT = Path(r); break\n\nassert DATA_ROOT is not None, \"RSNA 2022 dataset not found under /kaggle/input/**. Add it in 'Add data'.\"\n\nTRAIN_IMG_DIR = DATA_ROOT/\"train_images\"\nSEG_DIR = DATA_ROOT/\"segmentations\" # official folder name in competition data\nCSV_TRAIN = DATA_ROOT/\"train.csv\"\n\nprint(f\"[CHECK] DATA_ROOT: {DATA_ROOT}\")\nprint(f\"[CHECK] TRAIN_IMG_DIR exists: {TRAIN_IMG_DIR.exists()}\")\nprint(f\"[CHECK] SEG_DIR exists: {SEG_DIR.exists()}\")\nprint(f\"[CHECK] CSV_TRAIN exists: {CSV_TRAIN.exists()}\")\n\n# Enumerate studies (folders named by StudyInstanceUID)\nstudy_dirs = sorted([p for p in TRAIN_IMG_DIR.iterdir() if p.is_dir()])\nprint(f\"[EDA] #train studies found: {len(study_dirs)} (showing first 5)\")\nfor p in study_dirs[:5]: print(\" \", p.name)\n\n# Enumerate segmentation files (NIfTI .nii or .nii.gz)\nseg_files = sorted(list(SEG_DIR.glob(\"*.nii*\"))) if SEG_DIR.exists() else []\nprint(f\"[EDA] #segmentations found: {len(seg_files)} (showing first 5)\")\nfor f in seg_files[:5]: print(\" \", f.name)\n\n# Build mapping study_id -> segmentation path (by filename match)\nseg_map = {}\nfor f in seg_files:\n stem = f.name.split(\".nii\")[0]\n seg_map[stem] = f\ncoverage = sum(1 for p in study_dirs if p.name in seg_map)\nprint(f\"[EDA] Segmentation coverage: {coverage}/{len(study_dirs)} studies with masks\")\n\n# Load the label CSV (fracture labels; useful for cross-checks, not needed for stenosis)\nif CSV_TRAIN.exists():\n df_train = pd.read_csv(CSV_TRAIN)\n print(\"[EDA] train.csv shape:\", df_train.shape)\n display(df_train.head(3))\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T21:46:43.552815Z","iopub.execute_input":"2025-09-24T21:46:43.553374Z","iopub.status.idle":"2025-09-24T21:46:51.59299Z","shell.execute_reply.started":"2025-09-24T21:46:43.553347Z","shell.execute_reply":"2025-09-24T21:46:51.592327Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Cell 2: DICOM series loader, NIfTI loader, and geometry alignment ===\n\ndef pick_axial_series(study_folder: Path):\n \"\"\"\n Heuristic: prefer series path containing 'AX'/'AXIAL' and with the most DICOM files.\n Falls back to the subfolder with most files.\n \"\"\"\n cand = [p for p in study_folder.rglob(\"*\") if p.is_dir()]\n dcm_dirs = []\n for c in cand:\n try:\n n = sum(1 for _ in c.glob(\"*.dcm\"))\n if n >= 8: # skip tiny series\n hint = int((\"AX\" in c.name.upper()) or (\"AXIAL\" in c.name.upper()))\n dcm_dirs.append((hint, n, c))\n except Exception:\n pass\n if not dcm_dirs:\n return None\n dcm_dirs.sort(key=lambda t: (t[0], t[1]), reverse=True)\n return dcm_dirs[0][2]\n\ndef load_dicom_volume(series_dir: Path):\n reader = sitk.ImageSeriesReader()\n uids = reader.GetGDCMSeriesIDs(str(series_dir))\n if not uids: \n return None, None\n # choose UID with max slices\n best_uid, best_files = None, []\n for uid in uids:\n files = reader.GetGDCMSeriesFileNames(str(series_dir), uid)\n if len(files) > len(best_files):\n best_uid, best_files = uid, files\n reader.SetFileNames(best_files)\n img = reader.Execute()\n arr = sitk.GetArrayFromImage(img) # (Z,Y,X)\n return img, arr\n\ndef load_seg_nifti(seg_path: Path):\n seg_img = sitk.ReadImage(str(seg_path)) # supports .nii.gz\n seg_arr = sitk.GetArrayFromImage(seg_img) # (Z,Y,X)\n return seg_img, seg_arr\n\ndef resample_like(moving_img, reference_img, is_label=False):\n \"\"\"\n Resample 'moving_img' into 'reference_img' geometry (spacing, origin, direction).\n Nearest neighbor for labels, linear for images.\n \"\"\"\n interp = sitk.sitkNearestNeighbor if is_label else sitk.sitkLinear\n resamp = sitk.ResampleImageFilter()\n resamp.SetReferenceImage(reference_img)\n resamp.SetInterpolator(interp)\n resamp.SetTransform(sitk.Transform())\n resamp.SetDefaultPixelValue(0)\n out = resamp.Execute(moving_img)\n return out\n\ndef mid_sagittal_index(vol_arr):\n # pick mid X column\n return vol_arr.shape[2] // 2\n\ndef mid_axial_index(vol_arr):\n # pick mid Z slice\n return vol_arr.shape[0] // 2\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T21:47:07.170005Z","iopub.execute_input":"2025-09-24T21:47:07.170282Z","iopub.status.idle":"2025-09-24T21:47:07.181339Z","shell.execute_reply.started":"2025-09-24T21:47:07.170263Z","shell.execute_reply":"2025-09-24T21:47:07.180416Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Cell 2A: Robust series discovery (overrides) ===\nimport pydicom\nfrom collections import defaultdict\nimport math\n\ndef _series_score(first_file):\n \"\"\"Return (axial_hint, series_desc_hint) for scoring.\"\"\"\n axial_hint = 0\n desc_hint = 0\n try:\n ds = pydicom.dcmread(first_file, stop_before_pixels=True, specific_tags=[\n \"SeriesDescription\",\"ImageOrientationPatient\"\n ])\n sd = (ds.get(\"SeriesDescription\") or \"\").upper()\n if \"AX\" in sd or \"AXIAL\" in sd:\n axial_hint = 1\n # Orientation-based axial hint (normal ~ Z axis)\n iop = ds.get(\"ImageOrientationPatient\")\n if iop and len(iop)==6:\n # row, col direction cosines -> normal\n r = np.array(iop[:3], dtype=float); c = np.array(iop[3:], dtype=float)\n n = np.cross(r, c)\n if abs(n[2]) >= 0.8: # slice normal mostly along Z\n axial_hint = max(axial_hint, 1)\n # prefer typical CT image types\n if \"BONE\" in sd or \"STD\" in sd or \"SOFT\" in sd:\n desc_hint = 1\n except Exception:\n pass\n return axial_hint, desc_hint\n\ndef _gather_series(study_folder: Path, max_depth=2):\n \"\"\"Search study folder (and subfolders) for DICOM series; group by SeriesInstanceUID.\"\"\"\n dirs = [study_folder]\n if max_depth >= 1:\n dirs += [d for d in study_folder.glob(\"*\") if d.is_dir()]\n if max_depth >= 2:\n dirs += [d for d in study_folder.glob(\"*/*\") if d.is_dir()]\n series_map = defaultdict(list) # uid -> list of files\n first_file_for_uid = {}\n for d in dirs:\n # List .dcm files directly within this dir (no deeper to keep perf acceptable)\n dcm_files = sorted([str(x) for x in d.glob(\"*.dcm\")])\n if not dcm_files:\n continue\n # Group by SeriesInstanceUID\n for fp in dcm_files:\n try:\n ds = pydicom.dcmread(fp, stop_before_pixels=True, specific_tags=[\"SeriesInstanceUID\"])\n uid = str(ds.get(\"SeriesInstanceUID\"))\n if uid:\n series_map[uid].append(fp)\n if uid not in first_file_for_uid:\n first_file_for_uid[uid] = fp\n except Exception:\n continue\n # rank UIDs\n ranked = []\n for uid, files in series_map.items():\n axial_hint, desc_hint = _series_score(first_file_for_uid[uid])\n ranked.append((len(files), axial_hint, desc_hint, uid, sorted(files)))\n ranked.sort(reverse=True) # prefer more slices, axial, better desc\n return ranked # list of tuples\n\ndef load_best_series_any(study_folder: Path):\n \"\"\"\n Return (sitk_img, np_arr, files_used, uid) for the best series we can find.\n Tries direct SITK first; if none, uses pydicom grouping fallback.\n \"\"\"\n # Try SITK on study root\n uids = sitk.ImageSeriesReader.GetGDCMSeriesIDs(str(study_folder))\n best_tuple = None\n if uids:\n cand = []\n for uid in uids:\n files = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(str(study_folder), uid)\n # simple axial score by SeriesDescription (read first file)\n axial_hint, desc_hint = _series_score(files[0])\n cand.append((len(files), axial_hint, desc_hint, uid, files))\n cand.sort(reverse=True)\n best_tuple = cand[0]\n\n if best_tuple is None:\n # Fallback: scan subfolders, group by UID via pydicom\n ranked = _gather_series(study_folder, max_depth=2)\n if ranked:\n best_tuple = ranked[0]\n\n if best_tuple is None:\n return None, None, None, None\n\n nfiles, axial_hint, desc_hint, uid, files = best_tuple\n reader = sitk.ImageSeriesReader()\n reader.SetFileNames(files)\n img = reader.Execute()\n arr = sitk.GetArrayFromImage(img) # (Z,Y,X)\n print(f\"[LOAD] Series UID={uid} | slices={nfiles} | axial={axial_hint} | desc_hint={desc_hint}\")\n return img, arr, files, uid\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T21:51:48.370359Z","iopub.execute_input":"2025-09-24T21:51:48.370692Z","iopub.status.idle":"2025-09-24T21:51:49.119962Z","shell.execute_reply.started":"2025-09-24T21:51:48.370667Z","shell.execute_reply":"2025-09-24T21:51:49.119164Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Cell 3: Find a study with segmentation and visualize axial & sagittal overlays ===\n\ndef find_first_study_with_seg(study_dirs, seg_map):\n for p in study_dirs:\n if p.name in seg_map:\n return p, seg_map[p.name]\n return None, None\n\nstudy_dir, seg_path = find_first_study_with_seg(study_dirs, seg_map)\nassert study_dir is not None, \"No study with segmentation found. Check that 'segmentations' were added in 'Add data'.\"\n\nprint(f\"[PICK] study_id: {study_dir.name}\")\nseries_dir = pick_axial_series(study_dir)\nassert series_dir is not None, f\"No axial series found under {study_dir}\"\n\nprint(f\"[LOAD] DICOM series: {series_dir}\")\nct_img, ct_arr = load_dicom_volume(series_dir)\nassert ct_img is not None, \"Failed to load DICOM volume\"\n\nprint(f\"[LOAD] NIfTI seg: {seg_path.name}\")\nseg_img_mov, seg_arr_mov = load_seg_nifti(seg_path)\n\n# Resample segmentation into CT geometry if needed\nsame_geom = (ct_img.GetSize()==seg_img_mov.GetSize() and\n np.allclose(ct_img.GetSpacing(), seg_img_mov.GetSpacing(), atol=1e-3) and\n tuple(ct_img.GetDirection())==tuple(seg_img_mov.GetDirection()))\nif not same_geom:\n seg_img_ct = resample_like(seg_img_mov, ct_img, is_label=True)\n seg_arr = sitk.GetArrayFromImage(seg_img_ct)\n print(\"[INFO] Resampled segmentation to CT geometry.\")\nelse:\n seg_arr = seg_arr_mov\n\nprint(f\"[SHAPE] CT: {ct_arr.shape} | SEG: {seg_arr.shape} (Z,Y,X)\")\n\n# Choose display planes\nz_ax = mid_axial_index(ct_arr)\nx_sag = mid_sagittal_index(ct_arr)\n\n# Prepare slices\nax_ct = ct_arr[z_ax, :, :]\nax_lbl = seg_arr[z_ax, :, :]\nsag_ct = ct_arr[:, :, x_sag]\nsag_lbl = seg_arr[:, :, x_sag]\n\n# Window images and overlays\nax_vis = window_image(ax_ct, center=400, width=1800) # bone-friendly\nsag_vis = window_image(sag_ct, center=400, width=1800)\n\nax_ov = overlay_segmentation(ax_vis, ax_lbl)\nsag_ov = overlay_segmentation(sag_vis, sag_lbl)\n\n# Plot\nfig, axs = plt.subplots(1,2, figsize=(10,4.5))\naxs[0].imshow(ax_ov); axs[0].set_title(f\"Axial mid-slice (z={z_ax})\")\naxs[1].imshow(sag_ov); axs[1].set_title(f\"Sagittal mid-column (x={x_sag})\")\nfor a in axs: a.axis('off')\nplt.suptitle(f\"Study: {study_dir.name}\\nVertebra segmentation overlay (labels 1–7 = C1..C7)\")\nplt.show()\n\n# Quick label summary for this study\nvals, cnts = np.unique(seg_arr, return_counts=True)\npresent = {int(v): int(c) for v,c in zip(vals, cnts) if v>0}\nprint(\"[EDA] Labels present in segmentation (value: voxel_count) ->\", present)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T21:47:27.784203Z","iopub.execute_input":"2025-09-24T21:47:27.784513Z","iopub.status.idle":"2025-09-24T21:47:29.634149Z","shell.execute_reply.started":"2025-09-24T21:47:27.784492Z","shell.execute_reply":"2025-09-24T21:47:29.633028Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Cell 3A: Find a working CT+seg pair, resample seg, and show overlays ===\n\ndef find_working_pair_with_seg(study_dirs, seg_map, max_trials=50):\n trials = 0\n for p in study_dirs:\n if p.name not in seg_map: \n continue\n trials += 1\n if trials > max_trials:\n break\n try:\n ct_img, ct_arr, files, uid = load_best_series_any(p)\n if ct_img is None: \n print(f\"[SKIP] No series found for {p.name}\")\n continue\n seg_img_mov, seg_arr_mov = load_seg_nifti(seg_map[p.name])\n # resample seg to CT geometry if needed\n same_geom = (ct_img.GetSize()==seg_img_mov.GetSize() and\n np.allclose(ct_img.GetSpacing(), seg_img_mov.GetSpacing(), atol=1e-3) and\n tuple(ct_img.GetDirection())==tuple(seg_img_mov.GetDirection()))\n if not same_geom:\n seg_img_ct = resample_like(seg_img_mov, ct_img, is_label=True)\n seg_arr = sitk.GetArrayFromImage(seg_img_ct)\n print(f\"[INFO] Resampled seg -> CT geometry for {p.name}\")\n else:\n seg_arr = seg_arr_mov\n # sanity: must have some C1..C7 labels\n if not np.any(np.isin(seg_arr, np.arange(1,8))):\n print(f\"[SKIP] No cervical labels (1..7) in seg for {p.name}\")\n continue\n return p, ct_img, ct_arr, seg_arr\n except Exception as e:\n print(f\"[WARN] {p.name} failed: {e}\")\n continue\n return None, None, None, None\n\nstudy_dir2, ct_img2, ct_arr2, seg_arr2 = find_working_pair_with_seg(study_dirs, seg_map, max_trials=120)\nassert study_dir2 is not None, \"Could not find a study with both CT and cervical segmentation after 120 trials.\"\n\nprint(f\"[OK] Using study: {study_dir2.name} | CT shape={ct_arr2.shape} | SEG shape={seg_arr2.shape}\")\n\n# Choose display planes\nz_ax = mid_axial_index(ct_arr2)\nx_sag = mid_sagittal_index(ct_arr2)\n\n# Prepare slices\nax_ct = ct_arr2[z_ax, :, :]\nax_lbl = seg_arr2[z_ax, :, :]\nsag_ct = ct_arr2[:, :, x_sag]\nsag_lbl = seg_arr2[:, :, x_sag]\n\n# Window images and overlays\nax_vis = window_image(ax_ct, center=400, width=1800)\nsag_vis = window_image(sag_ct, center=400, width=1800)\n\nax_ov = overlay_segmentation(ax_vis, ax_lbl)\nsag_ov = overlay_segmentation(sag_vis, sag_lbl)\n\n# Plot\nfig, axs = plt.subplots(1,2, figsize=(10,4.5))\naxs[0].imshow(ax_ov); axs[0].set_title(f\"Axial mid-slice (z={z_ax})\")\naxs[1].imshow(sag_ov); axs[1].set_title(f\"Sagittal mid-column (x={x_sag})\")\nfor a in axs: a.axis('off')\nplt.suptitle(f\"Study: {study_dir2.name}\\nVertebra segmentation overlay (labels 1–7 = C1..C7)\")\nplt.show()\n\nvals, cnts = np.unique(seg_arr2, return_counts=True)\npresent = {int(v): int(c) for v,c in zip(vals, cnts) if v>0}\nprint(\"[EDA] Labels present (value: voxel_count) ->\", present)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T21:52:01.326207Z","iopub.execute_input":"2025-09-24T21:52:01.326526Z","iopub.status.idle":"2025-09-24T21:52:17.293181Z","shell.execute_reply.started":"2025-09-24T21:52:01.326503Z","shell.execute_reply":"2025-09-24T21:52:17.292088Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Recompute windows & overlays after the patch\nax_vis = window_image(ax_ct, center=400, width=1800)\nsag_vis = window_image(sag_ct, center=400, width=1800)\n\nax_ov = overlay_segmentation(ax_vis, ax_lbl)\nsag_ov = overlay_segmentation(sag_vis, sag_lbl)\n\nfig, axs = plt.subplots(1,2, figsize=(10,4.5))\naxs[0].imshow(ax_ov); axs[0].set_title(f\"Axial mid-slice (z={z_ax})\")\naxs[1].imshow(sag_ov); axs[1].set_title(f\"Sagittal mid-column (x={x_sag})\")\nfor a in axs: a.axis('off')\nplt.suptitle(f\"Study: {study_dir2.name}\\nVertebra segmentation overlay (labels 1–7=C1..C7; 8–19=T1..T12)\")\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T21:53:53.710484Z","iopub.execute_input":"2025-09-24T21:53:53.711294Z","iopub.status.idle":"2025-09-24T21:53:54.364295Z","shell.execute_reply.started":"2025-09-24T21:53:53.711271Z","shell.execute_reply":"2025-09-24T21:53:54.363336Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Cell 3P: Aspect-correct, orientation-aware display ===\n\ndef show_overlay(img2d, lbl2d, spacing_xy, title=\"\", flip_ud=False, flip_lr=False):\n \"\"\"\n img2d: 2D numpy (float in [0,1] recommended); lbl2d: int labels (same HW)\n spacing_xy: (sy, sx) in mm for the displayed plane\n flip_ud, flip_lr: visual-only flips to correct orientation\n \"\"\"\n v = img2d\n m = lbl2d\n if flip_ud: v, m = np.flipud(v), np.flipud(m)\n if flip_lr: v, m = np.fliplr(v), np.fliplr(m)\n\n ov = overlay_segmentation(v, m)\n sy, sx = float(spacing_xy[0]), float(spacing_xy[1])\n H, W = v.shape\n # extent maps pixel indices to physical mm so aspect is preserved\n extent = [0, W * sx, H * sy, 0] # x0, x1, y0, y1\n\n fig, ax = plt.subplots(figsize=(6, 6 * (H*sy)/(W*sx)))\n ax.imshow(ov, extent=extent, interpolation='nearest')\n ax.set_aspect('equal') # preserve mm aspect\n ax.set_title(title)\n ax.set_xlabel(\"mm\"); ax.set_ylabel(\"mm\")\n ax.axis('off')\n plt.tight_layout()\n plt.show()\n\n# Convenience wrappers for axial & sagittal using ct_img2 spacing:\n# ITK spacing order is (sx, sy, sz); array orders we show are:\n# axial: (Y,X) -> (sy, sx)\n# sagittal: (Z,Y) -> (sz, sy)\ndef show_axial(ax_img, ax_lbl, ct_img):\n sx, sy, sz = ct_img.GetSpacing()\n show_overlay(ax_img, ax_lbl, spacing_xy=(sy, sx),\n title=\"Axial mid-slice (aspect-correct)\")\n\ndef show_sagittal(sag_img, sag_lbl, ct_img, flip_upside_down=True):\n sx, sy, sz = ct_img.GetSpacing()\n # Many DICOMs render sagittal with superior at bottom; flip_ud=True fixes that visually\n show_overlay(sag_img, sag_lbl, spacing_xy=(sz, sy),\n title=\"Sagittal mid-column (aspect-correct)\",\n flip_ud=bool(flip_upside_down), flip_lr=False)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T21:55:28.205571Z","iopub.execute_input":"2025-09-24T21:55:28.206254Z","iopub.status.idle":"2025-09-24T21:55:28.215193Z","shell.execute_reply.started":"2025-09-24T21:55:28.206227Z","shell.execute_reply":"2025-09-24T21:55:28.214154Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Re-window just in case\nax_vis = window_image(ax_ct, center=400, width=1800)\nsag_vis = window_image(sag_ct, center=400, width=1800)\n\nshow_axial(ax_vis, ax_lbl, ct_img2)\nshow_sagittal(sag_vis, sag_lbl, ct_img2, flip_upside_down=True) # set False if your site already shows “upright”\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T21:55:31.249139Z","iopub.execute_input":"2025-09-24T21:55:31.250064Z","iopub.status.idle":"2025-09-24T21:55:31.830802Z","shell.execute_reply.started":"2025-09-24T21:55:31.250033Z","shell.execute_reply":"2025-09-24T21:55:31.829987Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Cell 4: Batch EDA — coverage stats + small overlay gallery ===\n\ndef study_has_cervical(seg_arr):\n return np.any(np.isin(seg_arr, np.arange(1,8))) # 1..7 = C1..C7\n\ndef overlay_panel(study_dir, seg_path, max_panels=4):\n try:\n series_dir = pick_axial_series(study_dir)\n ct_img, ct_arr = load_dicom_volume(series_dir)\n seg_img_mov, seg_arr_mov = load_seg_nifti(seg_path)\n # resample if needed\n same_geom = (ct_img.GetSize()==seg_img_mov.GetSize() and\n np.allclose(ct_img.GetSpacing(), seg_img_mov.GetSpacing(), atol=1e-3) and\n tuple(ct_img.GetDirection())==tuple(seg_img_mov.GetDirection()))\n seg_arr_ct = sitk.GetArrayFromImage(resample_like(seg_img_mov, ct_img, True)) if not same_geom else seg_arr_mov\n\n z = mid_axial_index(ct_arr); x = mid_sagittal_index(ct_arr)\n ax_vis = window_image(ct_arr[z]); ax_lbl = seg_arr_ct[z]\n sag_vis = window_image(ct_arr[:,:,x]); sag_lbl = seg_arr_ct[:,:,x]\n ax_ov = overlay_segmentation(ax_vis, ax_lbl)\n sag_ov = overlay_segmentation(sag_vis, sag_lbl)\n\n fig, axs = plt.subplots(1,2, figsize=(8.5,3.6))\n axs[0].imshow(ax_ov); axs[0].set_title(f\"Axial z={z}\")\n axs[1].imshow(sag_ov); axs[1].set_title(f\"Sagittal x={x}\")\n for a in axs: a.axis('off')\n plt.suptitle(f\"Study {study_dir.name} — Vertebra segmentation overlay\")\n plt.tight_layout(); plt.show()\n return True\n except Exception as e:\n print(f\"[WARN] Failed panel for {study_dir.name}: {e}\")\n return False\n\n# Coverage scan (first 300 studies for speed)\nhits = 0; cerv_hits = 0\nexamples = 0\nfor p in study_dirs[:300]:\n if p.name not in seg_map: continue\n seg_img, seg_arr = load_seg_nifti(seg_map[p.name])\n hits += 1\n if study_has_cervical(seg_arr): \n cerv_hits += 1\n if examples < 3:\n ok = overlay_panel(p, seg_map[p.name])\n if ok: examples += 1\n\nprint(f\"[EDA] In first 300 studies: {hits} have segmentation; {cerv_hits} contain C1–C7 labels.\")\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Cell 4A: Gallery (up to 4 studies) with crisp overlays ===\npicked = [study_dir2.name]\nshown = 0\nfor p in study_dirs:\n if p.name in picked or p.name not in seg_map:\n continue\n try:\n ct_img, ct_arr, files, uid = load_best_series_any(p)\n if ct_img is None: \n continue\n seg_img_mov, seg_arr_mov = load_seg_nifti(seg_map[p.name])\n same_geom = (ct_img.GetSize()==seg_img_mov.GetSize() and\n np.allclose(ct_img.GetSpacing(), seg_img_mov.GetSpacing(), atol=1e-3) and\n tuple(ct_img.GetDirection())==tuple(seg_img_mov.GetDirection()))\n seg_arr = sitk.GetArrayFromImage(resample_like(seg_img_mov, ct_img, True)) if not same_geom else seg_arr_mov\n if not np.any(np.isin(seg_arr, np.arange(1,8))):\n continue\n\n z = mid_axial_index(ct_arr); x = mid_sagittal_index(ct_arr)\n ax_vis = window_image(ct_arr[z]); ax_lbl = seg_arr[z]\n sag_vis = window_image(ct_arr[:,:,x]); sag_lbl = seg_arr[:,:,x]\n ax_ov = overlay_segmentation(ax_vis, ax_lbl)\n sag_ov = overlay_segmentation(sag_vis, sag_lbl)\n\n fig, axs = plt.subplots(1,2, figsize=(8.8,3.6))\n axs[0].imshow(ax_ov); axs[0].set_title(f\"Axial z={z}\")\n axs[1].imshow(sag_ov); axs[1].set_title(f\"Sagittal x={x}\")\n for a in axs: a.axis('off')\n plt.suptitle(f\"Study {p.name} — Vertebra segmentation overlay\")\n plt.tight_layout(); plt.show()\n\n picked.append(p.name); shown += 1\n if shown >= 3: break\n except Exception as e:\n print(f\"[WARN] Gallery skip {p.name}: {e}\")\n continue\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T21:58:26.254275Z","iopub.execute_input":"2025-09-24T21:58:26.255136Z","iopub.status.idle":"2025-09-24T21:58:51.724806Z","shell.execute_reply.started":"2025-09-24T21:58:26.255097Z","shell.execute_reply":"2025-09-24T21:58:51.723876Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Cell 5: Build per-study label presence table (C1..C7) ===\n\nrows = []\nfor p in study_dirs[:600]: # sample subset for speed; expand later\n sid = p.name\n has_seg = sid in seg_map\n row = {\"study_id\": sid, \"has_seg\": has_seg}\n if has_seg:\n try:\n seg_img, seg_arr = load_seg_nifti(seg_map[sid])\n for k in range(1,8):\n row[f\"C{k}\"] = int(np.any(seg_arr==k))\n except Exception as e:\n row.update({f\"C{k}\": -1 for k in range(1,8)})\n else:\n row.update({f\"C{k}\": 0 for k in range(1,8)})\n rows.append(row)\n\ndf_cov = pd.DataFrame(rows)\nprint(\"[EDA] per-study cervical label presence:\")\ndisplay(df_cov.head(10))\n\nprint(\"[EDA] counts:\\n\", df_cov[[f\"C{k}\" for k in range(1,8)]].sum())\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T21:59:01.839668Z","iopub.execute_input":"2025-09-24T21:59:01.839977Z","iopub.status.idle":"2025-09-24T21:59:36.520153Z","shell.execute_reply.started":"2025-09-24T21:59:01.839955Z","shell.execute_reply":"2025-09-24T21:59:36.519351Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Filter to “eligible” studies for measurement\neligible = []\nfor p in study_dirs:\n if p.name not in seg_map: \n continue\n seg_img, seg_arr = load_seg_nifti(seg_map[p.name])\n if np.any(np.isin(seg_arr, np.arange(1,8))): # any of C1..C7 present\n eligible.append(p)\nprint(f\"[FILTER] Eligible studies with cervical labels: {len(eligible)}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T22:01:32.415356Z","iopub.execute_input":"2025-09-24T22:01:32.415714Z","iopub.status.idle":"2025-09-24T22:03:32.747976Z","shell.execute_reply.started":"2025-09-24T22:01:32.41569Z","shell.execute_reply":"2025-09-24T22:03:32.745774Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def disc_mid_z_fallback(seg_arr, upper_id, lower_id):\n up = (seg_arr==upper_id); lo = (seg_arr==lower_id)\n if up.any() and lo.any():\n zu = int(np.median(np.argwhere(up)[:,0])); zl = int(np.median(np.argwhere(lo)[:,0]))\n return int(round((zu+zl)/2))\n # fallback: if only one present, use that vertebra's median Z\n if up.any():\n return int(np.median(np.argwhere(up)[:,0]))\n if lo.any():\n return int(np.median(np.argwhere(lo)[:,0]))\n return None\n\n# drop-in replacement inside Cell 5A:\ndef disc_mid_z(seg_arr, upper_id, lower_id):\n return disc_mid_z_fallback(seg_arr, upper_id, lower_id)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T22:03:43.791505Z","iopub.execute_input":"2025-09-24T22:03:43.79199Z","iopub.status.idle":"2025-09-24T22:03:43.800366Z","shell.execute_reply.started":"2025-09-24T22:03:43.791962Z","shell.execute_reply":"2025-09-24T22:03:43.799371Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Cell 5A: Per-level localization utilities ===\nCERV_IDS = np.arange(1, 8) # 1..7 = C1..C7 in RSNA2022\n\ndef level_presence(seg_arr):\n pres = {int(k): bool((seg_arr==k).any()) for k in CERV_IDS}\n return pres\n\ndef centroid_zy(mask3d):\n # returns (z,y) centroid in voxel coordinates for a 3D boolean mask\n idx = np.argwhere(mask3d)\n if idx.size == 0: return None\n z = int(np.median(idx[:,0])); y = int(np.median(idx[:,1]))\n return z, y\n\ndef disc_mid_z(seg_arr, upper_id, lower_id):\n # mid-z plane between two adjacent vertebra centroids (robust median)\n up = (seg_arr==upper_id)\n lo = (seg_arr==lower_id)\n cu = centroid_zy(up); cl = centroid_zy(lo)\n if cu is None or cl is None: return None\n return int(round((cu[0] + cl[0]) / 2))\n\ndef midsag_x_from_mask(seg_arr, lvl_id):\n # choose x column that traverses the vertebra mask centerline (robust median over occupied x)\n xs = np.argwhere(seg_arr==lvl_id)[:,2] if (seg_arr==lvl_id).any() else None\n if xs is None or xs.size==0: return seg_arr.shape[2]//2\n return int(np.median(xs))\n\n# Build a table of available C2..C7 levels with suggested planes\nrows = []\nfor lid in range(2,8): # C2..C7 planes defined by vertebra C{lid} and C{lid+1} (use prior if missing)\n rows.append(lid)\n\nperlevel = []\nfor lid in range(2,8): # define disc planes C2/3 ... C6/7\n upper, lower = lid, min(lid+1, 7)\n z_mid = disc_mid_z(seg_arr2, upper, lower) if lower!=upper else None\n x_mid = midsag_x_from_mask(seg_arr2, upper)\n perlevel.append({\"level\": f\"C{lid}/C{lower}\", \"upper\":upper, \"lower\":lower,\n \"z_mid\": z_mid, \"x_mid\": x_mid})\nperlevel_df = pd.DataFrame(perlevel)\nprint(\"[LOC] Proposed planes (nan if a vertebra label is missing):\")\ndisplay(perlevel_df)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T22:03:46.538814Z","iopub.execute_input":"2025-09-24T22:03:46.539128Z","iopub.status.idle":"2025-09-24T22:03:54.222471Z","shell.execute_reply.started":"2025-09-24T22:03:46.539105Z","shell.execute_reply":"2025-09-24T22:03:54.221582Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Cell 6B: Mask-guided APD & Torg (robust) ===\nBONE_TH = 200.0 # lower a bit for soft-tissue kernels; try 180–250 if needed\nZ_NEIGHBOR = 3 # median over ±3 slices around z_mid\n\ndef _posterior_body_from_mask(seg_band_zy, upper_id):\n \"\"\"\n seg_band_zy: (k, Y) int labels for a sagittal band (z-neighborhood at fixed x_mid)\n Return (y_vb_ant, y_vb_post) from the vertebral BODY mask (upper_id) projected along z.\n \"\"\"\n band = (seg_band_zy == int(upper_id)).astype(np.uint8) # (k,Y)\n proj = band.max(axis=0) # (Y,)\n ys = np.where(proj > 0)[0]\n if ys.size == 0:\n return None\n return int(ys.min()), int(ys.max())\n\ndef _lamina_anterior_y(ct_band_zy, y_start_post):\n \"\"\"\n Find anterior cortex of posterior elements: first contiguous bone cluster posterior to y_start_post.\n ct_band_zy: (k,Y) HU array; we binarize > BONE_TH and OR over k.\n \"\"\"\n bone = (ct_band_zy > BONE_TH)\n prof = bone.any(axis=0).astype(np.uint8) # (Y,)\n prof[:max(0, y_start_post+1)] = 0\n ys = np.where(prof > 0)[0]\n if ys.size == 0:\n return None\n return int(ys[0])\n\ndef measure_level_APD_Torg_mask_guided(ct_arr, ct_img, seg_arr, upper_id, z_mid, x_mid,\n z_neighborhood=Z_NEIGHBOR):\n \"\"\"\n Robust APD (mm) and Torg using vertebra mask guidance on a sagittal band.\n \"\"\"\n if z_mid is None or x_mid is None:\n return None\n sx, sy, sz = ct_img.GetSpacing() # (sx, sy, sz) -> use sy for AP distance in sagittal (Y)\n z0 = max(0, int(z_mid) - z_neighborhood)\n z1 = min(ct_arr.shape[0]-1, int(z_mid) + z_neighborhood)\n\n # build sagittal bands (k,Y)\n ct_band = ct_arr[z0:z1+1, :, int(x_mid)].astype(np.float32)\n seg_band = seg_arr[z0:z1+1, :, int(x_mid)].astype(np.int32)\n\n vb = _posterior_body_from_mask(seg_band, upper_id)\n if vb is None:\n return None\n y_vb_ant, y_vb_post = vb\n\n y_lam_ant = _lamina_anterior_y(ct_band, y_vb_post)\n if y_lam_ant is None or y_lam_ant <= y_vb_post:\n return None\n\n apd_mm = (y_lam_ant - y_vb_post) * float(sy)\n vb_ap_mm = max(0.0, (y_vb_post - y_vb_ant) * float(sy))\n torg = (apd_mm / vb_ap_mm) if vb_ap_mm > 0 else np.nan\n return {\"APD_mm\": float(apd_mm), \"Torg\": float(torg)}\n\ndef severity_from_apd_torg(apd_mm, torg):\n # Severe <10 mm; Moderate 10–13 mm; supportive Torg <0.8 (Pavlov/Torg). \n if apd_mm < 10.0: sev = 2\n elif apd_mm < 13.0: sev = 1\n else: sev = 0\n if not np.isnan(torg) and torg < 0.8 and 9.5 <= apd_mm <= 13.5:\n sev = min(2, sev+1)\n return sev\n\n# Recompute measurements for this study using mask-guided method\nmeas_rows = []\nfor r in perlevel:\n z_mid, x_mid, up = r[\"z_mid\"], r[\"x_mid\"], r[\"upper\"]\n if z_mid is None:\n meas_rows.append({**r, \"APD_mm\": np.nan, \"Torg\": np.nan, \"severity\": -1})\n continue\n m = measure_level_APD_Torg_mask_guided(ct_arr2, ct_img2, seg_arr2, up, z_mid, x_mid)\n if m is None:\n meas_rows.append({**r, \"APD_mm\": np.nan, \"Torg\": np.nan, \"severity\": -1})\n continue\n sev = severity_from_apd_torg(m[\"APD_mm\"], m[\"Torg\"])\n meas_rows.append({**r, **m, \"severity\": sev})\n\nmeas_df = pd.DataFrame(meas_rows)\nprint(\"[MEAS-v2] Per-level measurements (mask-guided):\")\ndisplay(meas_df)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T22:07:15.944908Z","iopub.execute_input":"2025-09-24T22:07:15.945209Z","iopub.status.idle":"2025-09-24T22:07:15.971383Z","shell.execute_reply.started":"2025-09-24T22:07:15.945188Z","shell.execute_reply":"2025-09-24T22:07:15.970252Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Cell 6C: Axial-plane mask-guided APD & Torg ===\nBONE_TH = 200.0 # soften if needed (180–250 works across kernels)\nX_BAND = 3 # median over x ∈ [x_mid±X_BAND]\nZ_NEIGH = 1 # small z jitter around disc plane if needed\n\ndef _axial_body_edges(ax_seg2d, vertebra_id):\n \"\"\"\n From axial segmentation at z=z_mid: for the given vertebra label,\n return (x_min, x_max, y_min, y_max) bounding box and the posterior\n body cortex y index (y_vb_post) estimated from the label mask.\n \"\"\"\n mask = (ax_seg2d == int(vertebra_id))\n if not mask.any():\n return None\n ys, xs = np.where(mask)\n y_min, y_max = ys.min(), ys.max()\n x_min, x_max = xs.min(), xs.max()\n # posterior body cortex approximated as the max Y within the vertebral body mask\n y_vb_post = int(y_max)\n return (x_min, x_max, y_min, y_max, y_vb_post)\n\ndef _axial_lamina_anterior_y(ax_ct2d, y_start_post, x_col):\n \"\"\"\n Scan posteriorly (in +Y) along the column x=x_col on axial CT,\n find the first bone cluster (> BONE_TH) posterior to y_start_post.\n Returns y_lam_ant or None.\n \"\"\"\n col = ax_ct2d[:, int(x_col)]\n bone = (col > BONE_TH).astype(np.uint8)\n bone[:max(0, y_start_post+1)] = 0 # zero anterior to/posterior body cortex\n ys = np.where(bone > 0)[0]\n return int(ys[0]) if ys.size else None\n\ndef measure_level_APD_Torg_axial(ct_arr, ct_img, seg_arr, up_id, z_mid, x_mid,\n x_band=X_BAND, z_neigh=Z_NEIGH):\n \"\"\"\n Compute APD (mm) and Torg on the **axial** slice near z_mid.\n AP = Y-axis in axial (array shape (Y,X)); lateral = X-axis.\n We use a small median over x ∈ [x_mid ± x_band] and z ∈ [z_mid ± z_neigh].\n \"\"\"\n if z_mid is None or x_mid is None:\n return None\n\n sx, sy, sz = ct_img.GetSpacing() # axial pixels: (sy, sx) → AP uses sy\n z0 = max(0, int(z_mid) - z_neigh)\n z1 = min(ct_arr.shape[0]-1, int(z_mid) + z_neigh)\n x0 = max(0, int(x_mid) - x_band)\n x1 = min(ct_arr.shape[2]-1, int(x_mid) + x_band)\n\n apd_vals, torg_vals = [], []\n\n for z in range(z0, z1+1):\n ax_ct = ct_arr[z].astype(np.float32) # (Y,X)\n ax_seg = seg_arr[z].astype(np.int32) # (Y,X)\n\n body_box = _axial_body_edges(ax_seg, up_id)\n if body_box is None:\n continue\n x_min, x_max, y_min, y_max, y_vb_post = body_box\n\n # VB anterior edge (rough): y_min from mask\n y_vb_ant = int(y_min)\n\n # scan across a small lateral band around x_mid\n y_lams = []\n for x in range(x0, x1+1):\n y_lam = _axial_lamina_anterior_y(ax_ct, y_vb_post, x)\n if y_lam is not None and y_lam > y_vb_post:\n y_lams.append(y_lam)\n\n if not y_lams:\n continue\n\n y_lam_ant = int(np.median(y_lams))\n apd_mm = max(0.0, (y_lam_ant - y_vb_post) * float(sy))\n vb_ap_mm = max(0.0, (y_vb_post - y_vb_ant) * float(sy))\n if vb_ap_mm <= 0:\n continue\n\n torg = apd_mm / vb_ap_mm\n apd_vals.append(apd_mm); torg_vals.append(torg)\n\n if not apd_vals:\n return None\n\n apd = float(np.median(apd_vals))\n torg = float(np.median(torg_vals)) if torg_vals else np.nan\n return {\"APD_mm\": apd, \"Torg\": torg}\n\ndef severity_from_apd_torg(apd_mm, torg):\n # Severe < 10 mm; Moderate 10–13 mm; Torg < 0.8 boosts near-boundaries\n if apd_mm < 10.0: sev = 2\n elif apd_mm < 13.0: sev = 1\n else: sev = 0\n if not np.isnan(torg) and torg < 0.8 and 9.5 <= apd_mm <= 13.5:\n sev = min(2, sev+1)\n return sev\n\n# Recompute measurements using axial method\nmeas_rows = []\nfor r in perlevel_df.to_dict(\"records\"):\n z_mid, x_mid, up = r[\"z_mid\"], r[\"x_mid\"], r[\"upper\"]\n if z_mid is None:\n meas_rows.append({**r, \"APD_mm\": np.nan, \"Torg\": np.nan, \"severity\": -1})\n continue\n m = measure_level_APD_Torg_axial(ct_arr2, ct_img2, seg_arr2, up, z_mid, x_mid)\n if m is None:\n meas_rows.append({**r, \"APD_mm\": np.nan, \"Torg\": np.nan, \"severity\": -1})\n continue\n sev = severity_from_apd_torg(m[\"APD_mm\"], m[\"Torg\"])\n meas_rows.append({**r, **m, \"severity\": sev})\n\nmeas_df = pd.DataFrame(meas_rows)\nprint(\"[MEAS-axial] Per-level measurements (mm) and severity:\")\ndisplay(meas_df)\n\n# === Patch: safe fallbacks for z_mid / x_mid when NaN ===\ndef safe_zmid(seg_arr, z_mid, upper_id, lower_id):\n # if z_mid is valid\n if z_mid == z_mid: # not NaN\n return int(round(z_mid))\n # fallback to median Z of whichever vertebra label exists\n for vid in (int(upper_id), int(lower_id)):\n m = (seg_arr == vid)\n if m.any():\n return int(np.median(np.argwhere(m)[:,0]))\n # last resort: center slice\n return seg_arr.shape[0] // 2\n\ndef safe_xmid(seg_arr, x_mid, upper_id):\n if x_mid == x_mid: # not NaN\n return int(round(x_mid))\n xs = np.argwhere(seg_arr == int(upper_id))[:,2] if (seg_arr == int(upper_id)).any() else None\n return int(np.median(xs)) if xs is not None and xs.size else seg_arr.shape[2] // 2\n\n# Re-run the measurement loop with safe fallbacks (uses measure_level_APD_Torg_axial from 6C)\nmeas_rows = []\nfor r in perlevel_df.to_dict(\"records\"):\n z_mid = safe_zmid(seg_arr2, r[\"z_mid\"], r[\"upper\"], r[\"lower\"])\n x_mid = safe_xmid(seg_arr2, r[\"x_mid\"], r[\"upper\"])\n up = int(r[\"upper\"])\n\n m = measure_level_APD_Torg_axial(ct_arr2, ct_img2, seg_arr2, up, z_mid, x_mid)\n if m is None:\n meas_rows.append({**r, \"z_mid\": z_mid, \"x_mid\": x_mid,\n \"APD_mm\": np.nan, \"Torg\": np.nan, \"severity\": -1})\n continue\n sev = severity_from_apd_torg(m[\"APD_mm\"], m[\"Torg\"])\n meas_rows.append({**r, \"z_mid\": z_mid, \"x_mid\": x_mid, **m, \"severity\": sev})\n\nmeas_df = pd.DataFrame(meas_rows)\nprint(\"[MEAS-axial] Per-level measurements with safe fallbacks:\")\ndisplay(meas_df)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T22:10:14.764894Z","iopub.execute_input":"2025-09-24T22:10:14.765248Z","iopub.status.idle":"2025-09-24T22:10:14.831831Z","shell.execute_reply.started":"2025-09-24T22:10:14.765225Z","shell.execute_reply":"2025-09-24T22:10:14.830704Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Cell 6D: Axial mask-guided APD/Torg with body-arch separation ===\n# Uses vertebra segmentation to isolate the vertebral BODY, then scans posteriorly to lamina/spinous cortex\n# on the union of all bony labels. Median across a small lateral band for robustness.\n\nimport numpy as np\nfrom scipy import ndimage as ndi\n\nBONE_TH = 180.0 # only used if you later want HU fallback; segmentation union is primary here\nX_BAND = 5 # median over x ∈ [x_mid ± X_BAND]\nZ_NEIGH = 1 # small z jitter around disc plane if needed\nOPEN_SZ = 5 # morphological opening size to remove thin posterior elements\n\ndef _axial_body_from_seg(ax_seg2d, vertebra_id):\n \"\"\"\n Return (y_vb_ant, y_vb_post, body_bbox) for the VERTEBRAL BODY only.\n Strategy: binary opening to delete thin arches; if still connected, keep anterior 65% of mask.\n \"\"\"\n mask = (ax_seg2d == int(vertebra_id))\n if not mask.any(): \n return None\n\n opened = ndi.binary_opening(mask, structure=np.ones((OPEN_SZ, OPEN_SZ)))\n body = opened.copy()\n if not body.any():\n # Fallback: cut posterior 35% of the original vertebra label\n ys = np.where(mask)[0]\n y_min, y_max = ys.min(), ys.max()\n y_cut = int(y_min + 0.65*(y_max - y_min))\n body = mask.copy()\n body[y_cut+1:, :] = 0\n\n ys_body, xs_body = np.where(body)\n if ys_body.size == 0:\n return None\n y_vb_ant = int(ys_body.min())\n y_vb_post = int(ys_body.max())\n x_min, x_max = int(xs_body.min()), int(xs_body.max())\n return y_vb_ant, y_vb_post, (x_min, x_max, int(ys_body.min()), int(ys_body.max()))\n\ndef _axial_lamina_anterior_y_from_seg(ax_seg2d, y_start_post, x_range):\n \"\"\"\n First posterior bony cortex behind the vertebral body:\n use the union of all labels (>0) so lamina/spinous processes are included.\n \"\"\"\n bone_union = (ax_seg2d > 0).astype(np.uint8)\n y_list = []\n for x in x_range:\n col = bone_union[:, int(x)].copy()\n col[:max(0, y_start_post+1)] = 0 # zero anterior to (and at) posterior body cortex\n ys = np.where(col > 0)[0]\n if ys.size:\n y_list.append(int(ys[0]))\n if not y_list:\n return None\n return int(np.median(y_list))\n\ndef measure_level_APD_Torg_axial(ct_arr, ct_img, seg_arr, up_id, z_mid, x_mid,\n x_band=X_BAND, z_neigh=Z_NEIGH):\n \"\"\"\n Compute APD (mm) and Torg on the axial plane near z_mid using segmentation only.\n AP is the Y-axis in axial; mm scaling uses sy from spacing.\n \"\"\"\n if z_mid is None or x_mid is None:\n return None\n\n sx, sy, sz = ct_img.GetSpacing() # axial px spacing: (sy, sx)\n z0 = max(0, int(z_mid) - z_neigh)\n z1 = min(ct_arr.shape[0]-1, int(z_mid) + z_neigh)\n x0 = max(0, int(x_mid) - x_band)\n x1 = min(ct_arr.shape[2]-1, int(x_mid) + x_band)\n\n apd_vals, torg_vals = [], []\n\n for z in range(z0, z1+1):\n ax_seg = seg_arr[z].astype(np.int32)\n\n vb = _axial_body_from_seg(ax_seg, up_id)\n if vb is None:\n continue\n y_vb_ant, y_vb_post, _bbox = vb\n\n # posterior bony cortex (lamina/spinous) from union of all labels\n y_lam_ant = _axial_lamina_anterior_y_from_seg(ax_seg, y_vb_post, range(x0, x1+1))\n if y_lam_ant is None or y_lam_ant <= y_vb_post:\n continue\n\n apd_mm = (y_lam_ant - y_vb_post) * float(sy)\n vb_ap_mm = max(0.0, (y_vb_post - y_vb_ant) * float(sy))\n if vb_ap_mm <= 0: \n continue\n\n torg = apd_mm / vb_ap_mm\n apd_vals.append(apd_mm)\n torg_vals.append(torg)\n\n if not apd_vals:\n return None\n apd = float(np.median(apd_vals))\n torg = float(np.median(torg_vals)) if torg_vals else np.nan\n return {\"APD_mm\": apd, \"Torg\": torg}\n\ndef severity_from_apd_torg(apd_mm, torg):\n # Severe <10 mm; Moderate 10–13 mm; Torg <0.8 supports up-binning near boundaries\n if apd_mm < 10.0: sev = 2\n elif apd_mm < 13.0: sev = 1\n else: sev = 0\n if not np.isnan(torg) and torg < 0.8 and 9.5 <= apd_mm <= 13.5:\n sev = min(2, sev+1)\n return sev\n\n# --- Recompute with safe fallbacks for z/x mids ---\ndef safe_zmid(seg_arr, z_mid, upper_id, lower_id):\n if z_mid == z_mid: # not NaN\n return int(round(z_mid))\n for vid in (int(upper_id), int(lower_id)):\n m = (seg_arr == vid)\n if m.any():\n return int(np.median(np.argwhere(m)[:,0]))\n return seg_arr.shape[0] // 2\n\ndef safe_xmid(seg_arr, x_mid, upper_id):\n if x_mid == x_mid:\n return int(round(x_mid))\n m = (seg_arr == int(upper_id))\n if m.any():\n xs = np.argwhere(m)[:,2]\n return int(np.median(xs))\n return seg_arr.shape[2] // 2\n\nmeas_rows = []\nfor r in perlevel_df.to_dict(\"records\"):\n z_mid = safe_zmid(seg_arr2, r[\"z_mid\"], r[\"upper\"], r[\"lower\"])\n x_mid = safe_xmid(seg_arr2, r[\"x_mid\"], r[\"upper\"])\n up = int(r[\"upper\"])\n m = measure_level_APD_Torg_axial(ct_arr2, ct_img2, seg_arr2, up, z_mid, x_mid)\n if m is None:\n meas_rows.append({**r, \"z_mid\": z_mid, \"x_mid\": x_mid, \"APD_mm\": np.nan, \"Torg\": np.nan, \"severity\": -1})\n continue\n sev = severity_from_apd_torg(m[\"APD_mm\"], m[\"Torg\"])\n meas_rows.append({**r, \"z_mid\": z_mid, \"x_mid\": x_mid, **m, \"severity\": sev})\n\nmeas_df = pd.DataFrame(meas_rows)\nprint(\"[MEAS-axial v2] Per-level measurements with body/arch separation:\")\ndisplay(meas_df)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T22:14:09.231141Z","iopub.execute_input":"2025-09-24T22:14:09.231995Z","iopub.status.idle":"2025-09-24T22:14:10.009287Z","shell.execute_reply.started":"2025-09-24T22:14:09.231968Z","shell.execute_reply":"2025-09-24T22:14:10.008397Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Patch: safe fallbacks for z_mid / x_mid when NaN ===\ndef safe_zmid(seg_arr, z_mid, upper_id, lower_id):\n # if z_mid is valid\n if z_mid == z_mid: # not NaN\n return int(round(z_mid))\n # fallback to median Z of whichever vertebra label exists\n for vid in (int(upper_id), int(lower_id)):\n m = (seg_arr == vid)\n if m.any():\n return int(np.median(np.argwhere(m)[:,0]))\n # last resort: center slice\n return seg_arr.shape[0] // 2\n\ndef safe_xmid(seg_arr, x_mid, upper_id):\n if x_mid == x_mid: # not NaN\n return int(round(x_mid))\n xs = np.argwhere(seg_arr == int(upper_id))[:,2] if (seg_arr == int(upper_id)).any() else None\n return int(np.median(xs)) if xs is not None and xs.size else seg_arr.shape[2] // 2\n\n# Re-run the measurement loop with safe fallbacks (uses measure_level_APD_Torg_axial from 6C)\nmeas_rows = []\nfor r in perlevel_df.to_dict(\"records\"):\n z_mid = safe_zmid(seg_arr2, r[\"z_mid\"], r[\"upper\"], r[\"lower\"])\n x_mid = safe_xmid(seg_arr2, r[\"x_mid\"], r[\"upper\"])\n up = int(r[\"upper\"])\n\n m = measure_level_APD_Torg_axial(ct_arr2, ct_img2, seg_arr2, up, z_mid, x_mid)\n if m is None:\n meas_rows.append({**r, \"z_mid\": z_mid, \"x_mid\": x_mid,\n \"APD_mm\": np.nan, \"Torg\": np.nan, \"severity\": -1})\n continue\n sev = severity_from_apd_torg(m[\"APD_mm\"], m[\"Torg\"])\n meas_rows.append({**r, \"z_mid\": z_mid, \"x_mid\": x_mid, **m, \"severity\": sev})\n\nmeas_df = pd.DataFrame(meas_rows)\nprint(\"[MEAS-axial] Per-level measurements with safe fallbacks:\")\ndisplay(meas_df)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T22:14:46.689416Z","iopub.execute_input":"2025-09-24T22:14:46.690285Z","iopub.status.idle":"2025-09-24T22:14:47.444423Z","shell.execute_reply.started":"2025-09-24T22:14:46.690249Z","shell.execute_reply":"2025-09-24T22:14:47.443277Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === QC (axial) — draw body posterior, lamina anterior, and APD on the chosen level ===\ndef qc_axial_level_v2(ct_arr, ct_img, seg_arr, level_row, x_band=X_BAND):\n z_mid, x_mid, up = int(level_row[\"z_mid\"]), int(level_row[\"x_mid\"]), int(level_row[\"upper\"])\n sx, sy, sz = ct_img.GetSpacing()\n ax_ct = ct_arr[z_mid].astype(np.float32)\n ax_seg = seg_arr[z_mid].astype(np.int32)\n\n vb = _axial_body_from_seg(ax_seg, up)\n if vb is None:\n print(\"No vertebral BODY on this slice.\")\n return\n y_vb_ant, y_vb_post, _ = vb\n\n x0 = max(0, x_mid - x_band); x1 = min(ax_ct.shape[1]-1, x_mid + x_band)\n y_lam_ant = _axial_lamina_anterior_y_from_seg(ax_seg, y_vb_post, range(x0, x1+1))\n\n img = window_image(ax_ct, 400, 1800)\n ov = overlay_segmentation(img, ax_seg)\n\n H,W = img.shape\n extent = [0, W*sx, H*sy, 0]\n fig, ax = plt.subplots(figsize=(6, 6*(H*sy)/(W*sx)))\n ax.imshow(ov, extent=extent); ax.set_aspect('equal'); ax.axis('off')\n\n # central column\n x_mm = x_mid * sx\n ax.axvline(x_mm, color=\"white\", lw=1.2, linestyle=\"--\")\n # posterior/anterior body, lamina\n ax.hlines([y_vb_ant*sy, y_vb_post*sy], x_mm-5, x_mm+5, colors=[\"cyan\",\"yellow\"], lw=2)\n if y_lam_ant is not None:\n ax.hlines([y_lam_ant*sy], x_mm-5, x_mm+5, colors=[\"red\"], lw=2)\n apd_mm = (y_lam_ant - y_vb_post) * sy\n ax.text(x_mm+6, (y_vb_post*sy + y_lam_ant*sy)/2,\n f\"APD ≈ {apd_mm:.1f} mm\", color=\"white\",\n bbox=dict(facecolor=\"black\", alpha=0.5), fontsize=10)\n\n ax.set_title(f\"Axial QC — {level_row['level']} (z={z_mid}, x={x_mid})\")\n plt.tight_layout(); plt.show()\n\nok = meas_df[meas_df[\"severity\"] >= 0]\nif len(ok):\n qc_axial_level_v2(ct_arr2, ct_img2, seg_arr2, ok.iloc[0])\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T22:14:14.981162Z","iopub.execute_input":"2025-09-24T22:14:14.98143Z","iopub.status.idle":"2025-09-24T22:14:15.404894Z","shell.execute_reply.started":"2025-09-24T22:14:14.981412Z","shell.execute_reply":"2025-09-24T22:14:15.403995Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === QC: draw axial APD column(s) for a measured level ===\ndef qc_axial_level(ct_arr, ct_img, seg_arr, level_row, x_band=X_BAND):\n z_mid, x_mid, up = int(level_row[\"z_mid\"]), int(level_row[\"x_mid\"]), int(level_row[\"upper\"])\n sx, sy, sz = ct_img.GetSpacing()\n ax_ct = ct_arr[z_mid].astype(np.float32)\n ax_seg = seg_arr[z_mid].astype(np.int32)\n\n body_box = _axial_body_edges(ax_seg, up)\n if body_box is None:\n print(\"No body mask on this slice.\")\n return\n x_min, x_max, y_min, y_max, y_vb_post = body_box\n y_vb_ant = int(y_min)\n\n x0 = max(0, x_mid - x_band); x1 = min(ax_ct.shape[1]-1, x_mid + x_band)\n\n # compute lamina positions used\n y_lams = []\n for x in range(x0, x1+1):\n y_lam = _axial_lamina_anterior_y(ax_ct, y_vb_post, x)\n if y_lam is not None and y_lam > y_vb_post:\n y_lams.append(y_lam)\n y_lam_ant = int(np.median(y_lams)) if y_lams else None\n\n img = window_image(ax_ct, 400, 1800)\n ov = overlay_segmentation(img, ax_seg)\n\n H,W = img.shape\n extent = [0, W*sx, H*sy, 0]\n fig, ax = plt.subplots(figsize=(6, 6*(H*sy)/(W*sx)))\n ax.imshow(ov, extent=extent); ax.set_aspect('equal'); ax.axis('off')\n\n # draw measurement column at x_mid\n x_mm = x_mid * sx\n ax.axvline(x_mm, color=\"white\", lw=1.5, linestyle=\"--\")\n\n # draw VB anterior/posterior and lamina lines along that column\n ax.hlines([y_vb_ant*sy, y_vb_post*sy], x_mm-4, x_mm+4, colors=[\"cyan\",\"yellow\"], lw=2)\n if y_lam_ant is not None:\n ax.hlines([y_lam_ant*sy], x_mm-4, x_mm+4, colors=[\"red\"], lw=2)\n\n apd_mm = (y_lam_ant - y_vb_post) * sy\n ax.text(x_mm+6, (y_vb_post*sy + y_lam_ant*sy)/2,\n f\"APD ≈ {apd_mm:.1f} mm\", color=\"white\",\n bbox=dict(facecolor=\"black\", alpha=0.5), fontsize=10)\n\n ax.set_title(f\"Axial QC — {level_row['level']} (z={z_mid}, x={x_mid})\")\n plt.tight_layout(); plt.show()\n\nok = meas_df[meas_df[\"severity\"] >= 0]\nif len(ok):\n qc_axial_level(ct_arr2, ct_img2, seg_arr2, ok.iloc[0])\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T22:12:28.917053Z","iopub.execute_input":"2025-09-24T22:12:28.917889Z","iopub.status.idle":"2025-09-24T22:12:29.382976Z","shell.execute_reply.started":"2025-09-24T22:12:28.917861Z","shell.execute_reply":"2025-09-24T22:12:29.381697Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Cell 6E: Axial APD/Torg with body-only mask + intensity lamina ===\nimport numpy as np\nfrom scipy import ndimage as ndi\n\nBONE_TH = 250.0 # try 220–300 depending on kernel\nX_BAND = 5 # lateral median band around x_mid\nZ_NEIGH = 1 # allow slight z jitter\nOPEN_SZ = 7 # opening kernel to delete thin arches\nKEEP_FRAC = 0.60 # keep anterior 60% of vertebra mask as \"body\"\nRUN_LEN = 3 # consecutive bone pixels to accept lamina\n\ndef _axial_body_only(ax_seg2d, vertebra_id):\n \"\"\"Return (y_vb_ant, y_vb_post) for vertebral BODY only.\"\"\"\n mask = (ax_seg2d == int(vertebra_id))\n if not mask.any(): \n return None\n # remove thin posterior elements\n opened = ndi.binary_opening(mask, structure=np.ones((OPEN_SZ, OPEN_SZ)))\n m = opened if opened.any() else mask.copy()\n ys, xs = np.where(m)\n y0, y1 = ys.min(), ys.max()\n # keep only anterior 60% to be body-only\n cut = int(y0 + KEEP_FRAC * (y1 - y0))\n body = m.copy(); body[cut+1:, :] = 0\n ys2 = np.where(body)[0]\n if ys2.size == 0: \n return None\n return int(ys2.min()), int(ys2.max())\n\ndef _lamina_y_from_intensity(ax_ct2d, y_start_post, x_range):\n \"\"\"Scan posteriorly for first robust bone run (intensity-based).\"\"\"\n H, W = ax_ct2d.shape\n y_min = max(0, y_start_post + 1)\n for y in range(y_min, H - RUN_LEN):\n # accept if ANY x in band has RUN_LEN consecutive bone pixels\n ok = False\n for x in x_range:\n col = ax_ct2d[y:y+RUN_LEN, int(x)]\n if np.all(col > BONE_TH):\n ok = True; break\n if ok:\n return y\n return None\n\ndef measure_level_APD_Torg_axial_v3(ct_arr, ct_img, seg_arr, up_id, z_mid, x_mid,\n x_band=X_BAND, z_neigh=Z_NEIGH):\n if z_mid is None or x_mid is None: \n return None\n sx, sy, sz = ct_img.GetSpacing()\n z0 = max(0, int(z_mid) - z_neigh)\n z1 = min(ct_arr.shape[0]-1, int(z_mid) + z_neigh)\n x0 = max(0, int(x_mid) - x_band)\n x1 = min(ct_arr.shape[2]-1, int(x_mid) + x_band)\n\n apd_vals, torg_vals = [], []\n for z in range(z0, z1+1):\n ax_ct = ct_arr[z].astype(np.float32)\n ax_seg = seg_arr[z].astype(np.int32)\n\n vb = _axial_body_only(ax_seg, up_id)\n if vb is None: \n continue\n y_vb_ant, y_vb_post = vb\n\n y_lam = _lamina_y_from_intensity(ax_ct, y_vb_post, range(x0, x1+1))\n if y_lam is None or y_lam <= y_vb_post:\n continue\n\n apd_mm = (y_lam - y_vb_post) * float(sy)\n vb_ap_mm = max(0.0, (y_vb_post - y_vb_ant) * float(sy))\n if vb_ap_mm <= 0:\n continue\n torg = apd_mm / vb_ap_mm\n apd_vals.append(apd_mm); torg_vals.append(torg)\n\n if not apd_vals:\n return None\n return {\"APD_mm\": float(np.median(apd_vals)),\n \"Torg\": float(np.median(torg_vals)) if torg_vals else np.nan}\n\ndef severity_from_apd_torg(apd_mm, torg):\n # Severe <10 mm; Moderate 10–13 mm; Torg<0.8 up-bins near the boundary.\n if apd_mm < 10.0: s = 2\n elif apd_mm < 13.0: s = 1\n else: s = 0\n if not np.isnan(torg) and torg < 0.8 and 9.5 <= apd_mm <= 13.5:\n s = min(2, s+1)\n return s\n\n# safe fallbacks (unchanged)\ndef safe_zmid(seg_arr, z_mid, upper_id, lower_id):\n if z_mid == z_mid: return int(round(z_mid))\n for vid in (int(upper_id), int(lower_id)):\n m = (seg_arr == vid)\n if m.any(): return int(np.median(np.argwhere(m)[:,0]))\n return seg_arr.shape[0] // 2\n\ndef safe_xmid(seg_arr, x_mid, upper_id):\n if x_mid == x_mid: return int(round(x_mid))\n m = (seg_arr == int(upper_id))\n if m.any():\n xs = np.argwhere(m)[:,2]\n return int(np.median(xs))\n return seg_arr.shape[2] // 2\n\n# --- recompute\nmeas_rows = []\nfor r in perlevel_df.to_dict(\"records\"):\n z_mid = safe_zmid(seg_arr2, r[\"z_mid\"], r[\"upper\"], r[\"lower\"])\n x_mid = safe_xmid(seg_arr2, r[\"x_mid\"], r[\"upper\"])\n up = int(r[\"upper\"])\n m = measure_level_APD_Torg_axial_v3(ct_arr2, ct_img2, seg_arr2, up, z_mid, x_mid)\n if m is None:\n meas_rows.append({**r, \"z_mid\": z_mid, \"x_mid\": x_mid,\n \"APD_mm\": np.nan, \"Torg\": np.nan, \"severity\": -1})\n continue\n sev = severity_from_apd_torg(m[\"APD_mm\"], m[\"Torg\"])\n meas_rows.append({**r, \"z_mid\": z_mid, \"x_mid\": x_mid, **m, \"severity\": sev})\nmeas_df = pd.DataFrame(meas_rows)\nprint(\"[MEAS-v3] Body-only + intensity lamina:\")\ndisplay(meas_df)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T22:17:56.472951Z","iopub.execute_input":"2025-09-24T22:17:56.47413Z","iopub.status.idle":"2025-09-24T22:17:57.349772Z","shell.execute_reply.started":"2025-09-24T22:17:56.474095Z","shell.execute_reply":"2025-09-24T22:17:57.34911Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# QC: axial overlay with v3 logic\ndef qc_axial_level_v3(ct_arr, ct_img, seg_arr, level_row, x_band=X_BAND):\n z_mid, x_mid, up = int(level_row[\"z_mid\"]), int(level_row[\"x_mid\"]), int(level_row[\"upper\"])\n sx, sy, sz = ct_img.GetSpacing()\n ax_ct = ct_arr[z_mid].astype(np.float32)\n ax_seg = seg_arr[z_mid].astype(np.int32)\n\n vb = _axial_body_only(ax_seg, up)\n if vb is None: \n print(\"No vertebral BODY on this slice.\"); return\n y_vb_ant, y_vb_post = vb\n x0 = max(0, x_mid - x_band); x1 = min(ax_ct.shape[1]-1, x_mid + x_band)\n y_lam = _lamina_y_from_intensity(ax_ct, y_vb_post, range(x0, x1+1))\n\n img = window_image(ax_ct, 400, 1800)\n ov = overlay_segmentation(img, ax_seg)\n H,W = img.shape\n extent = [0, W*sx, H*sy, 0]\n fig, ax = plt.subplots(figsize=(6, 6*(H*sy)/(W*sx)))\n ax.imshow(ov, extent=extent); ax.set_aspect('equal'); ax.axis('off')\n x_mm = x_mid * sx\n ax.axvline(x_mm, color=\"white\", lw=1.2, linestyle=\"--\")\n ax.hlines([y_vb_ant*sy, y_vb_post*sy], x_mm-6, x_mm+6, colors=[\"cyan\",\"yellow\"], lw=2)\n if y_lam is not None:\n ax.hlines([y_lam*sy], x_mm-6, x_mm+6, colors=[\"red\"], lw=2)\n apd_mm = (y_lam - y_vb_post) * sy\n ax.text(x_mm+6, (y_vb_post*sy + y_lam*sy)/2, f\"APD ≈ {apd_mm:.1f} mm\",\n color=\"white\", bbox=dict(facecolor=\"black\", alpha=0.5), fontsize=10)\n ax.set_title(f\"Axial QC — {level_row['level']} (z={z_mid}, x={x_mid})\")\n plt.tight_layout(); plt.show()\n\nok = meas_df[meas_df[\"severity\"] >= 0]\nif len(ok):\n qc_axial_level_v3(ct_arr2, ct_img2, seg_arr2, ok.iloc[0])\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T22:18:04.627329Z","iopub.execute_input":"2025-09-24T22:18:04.627611Z","iopub.status.idle":"2025-09-24T22:18:05.070099Z","shell.execute_reply.started":"2025-09-24T22:18:04.627591Z","shell.execute_reply":"2025-09-24T22:18:05.069146Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Cell 6F: Axial APD/Torg with ROI crop + artifact cleanup + robust stats ===\nimport numpy as np\nfrom scipy import ndimage as ndi\n\n# Tunables (you can tweak)\nBONE_TH = 240.0 # HU threshold for intensity fallback (220–300)\nOPEN_SZ = 7 # opening kernel to delete thin arches from vertebra label\nKEEP_FRAC = 0.65 # keep anterior 65% of vertebra label as \"body\"\nX_BAND = 5 # measure across x in [x_mid ± X_BAND]\nZ_NEIGH = 1 # z ∈ [z_mid ± Z_NEIGH]\nROI_Y_FRONT = 6 # mm anterior margin in ROI (beyond y_vb_ant)\nROI_Y_BACK = 25 # mm posterior margin in ROI (beyond y_vb_post)\nROI_X_MARGIN= 12 # mm lateral margin around x_mid\nAREA_MIN = 80 # px: remove small bone components (artifact spicules)\nROBUST_Q = 0.75 # take 75th percentile APD across band (robust to outliers)\nRUN_LEN = 3 # consecutive bone px to accept lamina\n\ndef _axial_body_only(ax_seg2d, vertebra_id):\n m = (ax_seg2d == int(vertebra_id))\n if not m.any(): return None\n opened = ndi.binary_opening(m, structure=np.ones((OPEN_SZ, OPEN_SZ)))\n base = opened if opened.any() else m\n ys, xs = np.where(base)\n y0, y1 = ys.min(), ys.max()\n cut = int(y0 + KEEP_FRAC*(y1 - y0))\n body = base.copy(); body[cut+1:, :] = 0\n ys2 = np.where(body)[0]\n if ys2.size == 0: return None\n return int(ys2.min()), int(ys2.max())\n\ndef _px_from_mm(mm, spacing): # helper\n return int(round(mm / float(spacing)))\n\ndef _clean_bone_union(ax_seg2d):\n # union all labels >0; remove tiny components\n bone = (ax_seg2d > 0).astype(np.uint8)\n lab, num = ndi.label(bone)\n if num == 0: return bone\n sizes = np.bincount(lab.ravel())\n kill = np.where(sizes < AREA_MIN)[0]\n bone[np.isin(lab, kill)] = 0\n return bone.astype(np.uint8)\n\ndef _lamina_y_from_union(ax_bone2d, y_start_post, x_range):\n # first robust bone posterior to body, requiring RUN_LEN consecutive bone px\n H, W = ax_bone2d.shape\n y_min = max(0, y_start_post + 1)\n for y in range(y_min, H-RUN_LEN):\n ok = False\n for x in x_range:\n if np.all(ax_bone2d[y:y+RUN_LEN, int(x)] > 0):\n ok = True; break\n if ok:\n return y\n return None\n\ndef measure_level_APD_Torg_axial_v4(ct_arr, ct_img, seg_arr, up_id, z_mid, x_mid):\n if z_mid is None or x_mid is None: return None\n sx, sy, sz = ct_img.GetSpacing() # axial pixels: (sy, sx)\n z0 = max(0, int(z_mid) - Z_NEIGH)\n z1 = min(ct_arr.shape[0]-1, int(z_mid) + Z_NEIGH)\n\n # mm → px ROI margins\n dy_front = _px_from_mm(ROI_Y_FRONT, sy)\n dy_back = _px_from_mm(ROI_Y_BACK, sy)\n dxm = _px_from_mm(ROI_X_MARGIN, sx)\n\n apd_all, torg_all = [], []\n\n for z in range(z0, z1+1):\n ax_ct = ct_arr[z].astype(np.float32)\n ax_seg = seg_arr[z].astype(np.int32)\n\n # 1) BODY-only from label\n vb = _axial_body_only(ax_seg, up_id)\n if vb is None: \n continue\n y_vb_ant, y_vb_post = vb\n\n # 2) ROI crop (“cone in”)\n H, W = ax_ct.shape\n y0 = max(0, y_vb_ant - dy_front)\n y1 = min(H-1, y_vb_post + dy_back)\n x0 = max(0, int(x_mid) - dxm)\n x1 = min(W-1, int(x_mid) + dxm)\n\n ax_ct_roi = ax_ct[y0:y1+1, x0:x1+1]\n ax_seg_roi = ax_seg[y0:y1+1, x0:x1+1]\n y_vb_ant_r = y_vb_ant - y0\n y_vb_post_r= y_vb_post - y0\n x_mid_r = int(x_mid) - x0\n\n # 3) bone union cleanup in ROI\n bone_u = _clean_bone_union(ax_seg_roi) # 0/1\n\n # 4) measure across lateral band (robust to one bad column)\n xb0 = max(0, x_mid_r - X_BAND)\n xb1 = min(ax_ct_roi.shape[1]-1, x_mid_r + X_BAND)\n apds, torgs = [], []\n\n for x in range(xb0, xb1+1):\n y_lam = _lamina_y_from_union(bone_u, y_vb_post_r, [x])\n if y_lam is None or y_lam <= y_vb_post_r: \n continue\n apd_mm = (y_lam - y_vb_post_r) * float(sy)\n vb_ap_mm = max(0.0, (y_vb_post_r - y_vb_ant_r) * float(sy))\n if vb_ap_mm <= 0: \n continue\n apds.append(apd_mm); torgs.append(apd_mm / vb_ap_mm)\n\n if len(apds) == 0:\n continue\n\n apd_rob = float(np.quantile(apds, ROBUST_Q))\n torg_rob= float(np.quantile(torgs, ROBUST_Q)) if len(torgs) else np.nan\n apd_all.append(apd_rob); torg_all.append(torg_rob)\n\n if len(apd_all) == 0:\n return None\n apd = float(np.median(apd_all))\n torg = float(np.median(torg_all)) if len(torg_all) else np.nan\n return {\"APD_mm\": apd, \"Torg\": torg}\n\ndef severity_from_apd_torg(apd_mm, torg):\n # Severe <10 mm; Moderate 10–13 mm; Torg<0.8 supports up-binning near boundary\n if apd_mm < 10.0: s = 2\n elif apd_mm < 13.0: s = 1\n else: s = 0\n if not np.isnan(torg) and torg < 0.8 and 9.5 <= apd_mm <= 13.5:\n s = min(2, s+1)\n return s\n\n# Safe fallbacks (same as before)\ndef safe_zmid(seg_arr, z_mid, upper_id, lower_id):\n if z_mid == z_mid: return int(round(z_mid))\n for vid in (int(upper_id), int(lower_id)):\n m = (seg_arr == vid)\n if m.any(): return int(np.median(np.argwhere(m)[:,0]))\n return seg_arr.shape[0] // 2\n\ndef safe_xmid(seg_arr, x_mid, upper_id):\n if x_mid == x_mid: return int(round(x_mid))\n m = (seg_arr == int(upper_id))\n if m.any():\n xs = np.argwhere(m)[:,2]\n return int(np.median(xs))\n return seg_arr.shape[2] // 2\n\n# --- recompute for this study with v4 ---\nmeas_rows = []\nfor r in perlevel_df.to_dict(\"records\"):\n z_mid = safe_zmid(seg_arr2, r[\"z_mid\"], r[\"upper\"], r[\"lower\"])\n x_mid = safe_xmid(seg_arr2, r[\"x_mid\"], r[\"upper\"])\n up = int(r[\"upper\"])\n m = measure_level_APD_Torg_axial_v4(ct_arr2, ct_img2, seg_arr2, up, z_mid, x_mid)\n if m is None:\n meas_rows.append({**r, \"z_mid\": z_mid, \"x_mid\": x_mid, \"APD_mm\": np.nan, \"Torg\": np.nan, \"severity\": -1})\n continue\n sev = severity_from_apd_torg(m[\"APD_mm\"], m[\"Torg\"])\n meas_rows.append({**r, \"z_mid\": z_mid, \"x_mid\": x_mid, **m, \"severity\": sev})\n\nmeas_df = pd.DataFrame(meas_rows)\nprint(\"[MEAS-v4] ROI + cleanup + robust APD:\")\ndisplay(meas_df)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T22:22:11.561968Z","iopub.execute_input":"2025-09-24T22:22:11.562772Z","iopub.status.idle":"2025-09-24T22:22:12.552666Z","shell.execute_reply.started":"2025-09-24T22:22:11.562744Z","shell.execute_reply":"2025-09-24T22:22:12.551967Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === QC for v4: compact, aspect-correct panel ===\ndef qc_axial_level_v4(ct_arr, ct_img, seg_arr, level_row):\n z_mid, x_mid, up = int(level_row[\"z_mid\"]), int(level_row[\"x_mid\"]), int(level_row[\"upper\"])\n sx, sy, sz = ct_img.GetSpacing()\n ax_ct = ct_arr[z_mid].astype(np.float32)\n ax_seg = seg_arr[z_mid].astype(np.int32)\n\n vb = _axial_body_only(ax_seg, up)\n if vb is None: \n print(\"No vertebral BODY on this slice.\"); return\n y_vb_ant, y_vb_post = vb\n\n # ROI as above\n dy_front = _px_from_mm(ROI_Y_FRONT, sy); dy_back = _px_from_mm(ROI_Y_BACK, sy); dxm = _px_from_mm(ROI_X_MARGIN, sx)\n H,W = ax_ct.shape\n y0 = max(0, y_vb_ant - dy_front); y1 = min(H-1, y_vb_post + dy_back)\n x0 = max(0, int(x_mid) - dxm); x1 = min(W-1, int(x_mid) + dxm)\n ct_roi = ax_ct[y0:y1+1, x0:x1+1]\n seg_roi = ax_seg[y0:y1+1, x0:x1+1]\n bone_u = _clean_bone_union(seg_roi)\n\n img = window_image(ct_roi, center=500, width=2200) # boney display\n ov = overlay_segmentation(img, seg_roi)\n\n # measure at center column after cleanup (for display)\n y_vb_post_r = y_vb_post - y0\n x_mid_r = int(x_mid) - x0\n y_lam = _lamina_y_from_union(bone_u, y_vb_post_r, [x_mid_r])\n\n H2,W2 = img.shape\n extent = [0, W2*sx, H2*sy, 0]\n fig, ax = plt.subplots(figsize=(5.8, 5.8*(H2*sy)/(W2*sx)))\n ax.imshow(ov, extent=extent); ax.set_aspect('equal'); ax.axis('off')\n ax.axvline(x_mid_r*sx, color=\"white\", lw=1.2, linestyle=\"--\")\n ax.hlines([ (y_vb_post_r)*sy ], x_mid_r*sx-6, x_mid_r*sx+6, colors=\"yellow\", lw=2, label=\"VB posterior\")\n if y_lam is not None:\n ax.hlines([ y_lam*sy ], x_mid_r*sx-6, x_mid_r*sx+6, colors=\"red\", lw=2, label=\"Lamina anterior\")\n apd_mm = (y_lam - y_vb_post_r) * sy\n ax.text(x_mid_r*sx+6, (y_vb_post_r*sy + y_lam*sy)/2,\n f\"APD ≈ {apd_mm:.1f} mm\", color=\"white\",\n bbox=dict(facecolor=\"black\", alpha=0.55), fontsize=10)\n ax.set_title(f\"Axial QC — {level_row['level']} (coned-in)\")\n plt.tight_layout(); plt.show()\n\nok = meas_df[meas_df[\"severity\"] >= 0]\nif len(ok):\n qc_axial_level_v4(ct_arr2, ct_img2, seg_arr2, ok.iloc[0])\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T22:22:31.469127Z","iopub.execute_input":"2025-09-24T22:22:31.46996Z","iopub.status.idle":"2025-09-24T22:22:31.837865Z","shell.execute_reply.started":"2025-09-24T22:22:31.469935Z","shell.execute_reply":"2025-09-24T22:22:31.836959Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Cell 6G: Midline-refined axial APD/Torg + dual-view QC ===\nimport numpy as np\nfrom scipy import ndimage as ndi\n\n# Tunables\nBONE_TH = 240.0\nOPEN_SZ = 7\nKEEP_FRAC = 0.65 # keep anterior 65% of vertebra mask as \"body-only\"\nX_BAND = 5 # band around the midline for robust stats\nZ_NEIGH = 1 # z jitter around disc mid-plane\nROI_Y_FRONT = 10 # mm anterior margin in ROI (beyond VB anterior)\nROI_Y_BACK = 30 # mm posterior margin in ROI (beyond VB posterior)\nROI_X_MARG = 20 # mm lateral margin around x_mid\nAREA_MIN = 120 # px: remove tiny bone components in ROI\nRUN_LEN = 3 # consecutive bone px (de-streaking)\nROBUST_Q = 0.75 # percentile across band (robust APD)\nMED_FILT = (3, 1) # median filter kernel for ROI (reduce streaks)\n\ndef _px_from_mm(mm, sp): return int(round(mm / float(sp)))\n\ndef _body_only_mask(ax_seg2d, vid):\n m = (ax_seg2d == int(vid))\n if not m.any(): return None\n opened = ndi.binary_opening(m, structure=np.ones((OPEN_SZ, OPEN_SZ)))\n base = opened if opened.any() else m\n ys, xs = np.where(base)\n y0, y1 = ys.min(), ys.max()\n cut = int(y0 + KEEP_FRAC * (y1 - y0))\n body = base.copy()\n body[cut+1:, :] = 0\n return body if body.any() else None\n\ndef _refine_x_mid_from_bodies(ax_seg2d, up_id, lo_id, x_default):\n xs = []\n for vid in (int(up_id), int(lo_id)):\n m = (ax_seg2d == vid)\n if m.any():\n xs.append(int(np.median(np.where(m)[1])))\n return int(np.mean(xs)) if xs else int(x_default)\n\ndef _clean_bone_union(ax_seg2d):\n bone = (ax_seg2d > 0).astype(np.uint8)\n lab, num = ndi.label(bone)\n if num == 0: return bone\n sizes = np.bincount(lab.ravel())\n kill = np.where(sizes < AREA_MIN)[0]\n bone[np.isin(lab, kill)] = 0\n return bone.astype(np.uint8)\n\ndef _first_bone_run_y(bone_bin2d, y_from, x):\n H, W = bone_bin2d.shape\n y0 = max(0, int(y_from) + 1)\n for y in range(y0, H - RUN_LEN):\n if np.all(bone_bin2d[y:y+RUN_LEN, int(x)] > 0):\n return y\n return None\n\ndef _safe_zmid(seg_arr, z_mid, up_id, lo_id):\n if z_mid == z_mid: return int(round(z_mid))\n for vid in (int(up_id), int(lo_id)):\n m = (seg_arr == vid)\n if m.any(): return int(np.median(np.argwhere(m)[:,0]))\n return seg_arr.shape[0] // 2\n\ndef _safe_xmid(seg_arr, x_mid, up_id, lo_id, z_mid):\n if x_mid == x_mid: return int(round(x_mid))\n # try centroid at z_mid from both vertebrae; else default to image center\n if 0 <= z_mid < seg_arr.shape[0]:\n return _refine_x_mid_from_bodies(seg_arr[z_mid], up_id, lo_id, seg_arr.shape[2]//2)\n return seg_arr.shape[2] // 2\n\ndef measure_axial_midline(ct_arr, ct_img, seg_arr, up_id, lo_id, z_mid, x_mid):\n \"\"\"Returns dict with APD_mid (mm), APD_band (mm), Torg_mid, Torg_band, and plotting info.\"\"\"\n sx, sy, sz = ct_img.GetSpacing()\n z_mid = _safe_zmid(seg_arr, z_mid, up_id, lo_id)\n x_mid = _safe_xmid(seg_arr, x_mid, up_id, lo_id, z_mid)\n\n z0 = max(0, z_mid - Z_NEIGH)\n z1 = min(ct_arr.shape[0]-1, z_mid + Z_NEIGH)\n\n # compute on each nearby z, collect stats\n apd_mid_list, torg_mid_list, apd_band_list, torg_band_list = [], [], [], []\n # store last ROI slices for QC\n qc_payload = None\n\n for z in range(z0, z1+1):\n ax_ct = ct_arr[z].astype(np.float32)\n ax_seg = seg_arr[z].astype(np.int32)\n\n body_mask = _body_only_mask(ax_seg, up_id)\n if body_mask is None: \n continue\n\n ys, xs = np.where(body_mask)\n y_vb_ant = int(ys.min()); y_vb_post = int(ys.max())\n\n # ROI (cone-in, but wide enough to include entire spine)\n dy_f = _px_from_mm(ROI_Y_FRONT, sy)\n dy_b = _px_from_mm(ROI_Y_BACK, sy)\n dxm = _px_from_mm(ROI_X_MARG, sx)\n\n H,W = ax_ct.shape\n y0 = max(0, y_vb_ant - dy_f); y1 = min(H-1, y_vb_post + dy_b)\n x0 = max(0, x_mid - dxm); x1 = min(W-1, x_mid + dxm)\n\n ct_roi = ax_ct[y0:y1+1, x0:x1+1]\n seg_roi = ax_seg[y0:y1+1, x0:x1+1]\n body_roi= body_mask[y0:y1+1, x0:x1+1]\n bone_u = _clean_bone_union(seg_roi)\n\n # slight de-streak\n ct_roi_f = ndi.median_filter(ct_roi, size=MED_FILT)\n\n # midline & band in ROI coords\n x_mid_r = int(x_mid) - x0\n y_vb_ant_r, y_vb_post_r = y_vb_ant - y0, y_vb_post - y0\n xb0 = max(0, x_mid_r - X_BAND); xb1 = min(ct_roi.shape[1]-1, x_mid_r + X_BAND)\n\n # --- MIDLINE measurement\n y_lam_mid = _first_bone_run_y(bone_u, y_vb_post_r, x_mid_r)\n if y_lam_mid is not None and y_lam_mid > y_vb_post_r:\n apd_mid = (y_lam_mid - y_vb_post_r) * float(sy)\n vb_ap = max(0.0, (y_vb_post_r - y_vb_ant_r) * float(sy))\n if vb_ap > 0:\n apd_mid_list.append(apd_mid)\n torg_mid_list.append(apd_mid / vb_ap)\n\n # --- BAND-ROBUST measurement\n apds, torgs = [], []\n for x in range(xb0, xb1+1):\n y_lam = _first_bone_run_y(bone_u, y_vb_post_r, x)\n if y_lam is None or y_lam <= y_vb_post_r: \n continue\n apd = (y_lam - y_vb_post_r) * float(sy)\n vb_ap = max(0.0, (y_vb_post_r - y_vb_ant_r) * float(sy))\n if vb_ap > 0:\n apds.append(apd); torgs.append(apd / vb_ap)\n if apds:\n apd_band_list.append(float(np.quantile(apds, ROBUST_Q)))\n torg_band_list.append(float(np.quantile(torgs, ROBUST_Q)))\n\n qc_payload = dict(z=z, y0=y0, y1=y1, x0=x0, x1=x1, \n x_mid_r=x_mid_r, y_vb_ant_r=y_vb_ant_r, y_vb_post_r=y_vb_post_r,\n y_lam_mid=y_lam_mid, ct_roi=ct_roi_f, seg_roi=seg_roi)\n\n if not (apd_mid_list or apd_band_list):\n return None\n\n res = dict(\n APD_mid_mm = np.median(apd_mid_list) if apd_mid_list else np.nan,\n APD_band_mm= np.median(apd_band_list) if apd_band_list else np.nan,\n Torg_mid = np.median(torg_mid_list) if torg_mid_list else np.nan,\n Torg_band = np.median(torg_band_list) if torg_band_list else np.nan,\n z_mid_used = z_mid,\n x_mid_used = x_mid,\n qc=qc_payload\n )\n return res\n\ndef severity_from_apd(apd_mm):\n # severe <10, moderate 10–13, else none/mild\n if np.isnan(apd_mm): return -1\n return 2 if apd_mm < 10.0 else (1 if apd_mm < 13.0 else 0)\n\n# ---- run over this study\nrows = []\nfor r in perlevel_df.to_dict(\"records\"):\n m = measure_axial_midline(ct_arr2, ct_img2, seg_arr2, r[\"upper\"], r[\"lower\"], r[\"z_mid\"], r[\"x_mid\"])\n if m is None:\n rows.append({**r, \"APD_mid_mm\": np.nan, \"APD_band_mm\": np.nan, \"Torg_mid\": np.nan, \"Torg_band\": np.nan,\n \"severity_mid\": -1, \"severity_band\": -1})\n else:\n rows.append({**r, **{k:v for k,v in m.items() if k!=\"qc\"},\n \"severity_mid\": severity_from_apd(m[\"APD_mid_mm\"]),\n \"severity_band\": severity_from_apd(m[\"APD_band_mm\"])})\nmeas_df = pd.DataFrame(rows)\nprint(\"[MEAS-midline] APD (midline & band) and severities:\")\ndisplay(meas_df)\n\n# ---- Dual-view QC (full + ROI) for the first measurable level\ndef qc_dual_full_and_roi(ct_arr, ct_img, seg_arr, level_row, qc_payload,\n window_full=(400,1800), window_roi=(500,2200)):\n sx, sy, sz = ct_img.GetSpacing()\n z = qc_payload[\"z\"]; y0,y1,x0,x1 = qc_payload[\"y0\"],qc_payload[\"y1\"],qc_payload[\"x0\"],qc_payload[\"x1\"]\n x_mid_r = qc_payload[\"x_mid_r\"]; y_vb_ant_r=qc_payload[\"y_vb_ant_r\"]; y_vb_post_r=qc_payload[\"y_vb_post_r\"]\n y_lam_mid = qc_payload[\"y_lam_mid\"]; ct_roi = qc_payload[\"ct_roi\"]; seg_roi = qc_payload[\"seg_roi\"]\n\n # full-FOV left\n full = window_image(ct_arr[z], *window_full); seg_full = seg_arr[z]\n fig, axs = plt.subplots(1,2, figsize=(11,5.5))\n\n extent_full = [0, full.shape[1]*sx, full.shape[0]*sy, 0]\n axs[0].imshow(overlay_segmentation(full, seg_full), extent=extent_full); axs[0].set_aspect('equal'); axs[0].axis('off')\n axs[0].axvline(meas_df.loc[0,\"x_mid\"], color=\"white\", lw=1.2, linestyle=\"--\")\n axs[0].set_title(f\"Full axial (z={z})\")\n\n # ROI right\n roi = window_image(ct_roi, *window_roi)\n extent_roi = [0, roi.shape[1]*sx, roi.shape[0]*sy, 0]\n axs[1].imshow(overlay_segmentation(roi, seg_roi), extent=extent_roi); axs[1].set_aspect('equal'); axs[1].axis('off')\n x_mm = x_mid_r * sx\n axs[1].axvline(x_mm, color=\"white\", lw=1.2, linestyle=\"--\")\n axs[1].hlines([y_vb_ant_r*sy, y_vb_post_r*sy], x_mm-8, x_mm+8, colors=[\"cyan\",\"yellow\"], lw=2)\n if y_lam_mid is not None:\n axs[1].hlines([y_lam_mid*sy], x_mm-8, x_mm+8, colors=[\"red\"], lw=2)\n apd = (y_lam_mid - y_vb_post_r) * sy\n axs[1].text(x_mm+6, (y_vb_post_r*sy + y_lam_mid*sy)/2,\n f\"APD(mid) ≈ {apd:.1f} mm\", color=\"white\",\n bbox=dict(facecolor=\"black\", alpha=0.55), fontsize=10)\n axs[1].set_title(f\"Coned-in ROI\")\n\n plt.suptitle(f\"QC — {level_row['level']} (midline & ROI)\", y=0.98, fontsize=12)\n plt.tight_layout(); plt.show()\n\n# pick first measurable level\nidx = meas_df.index[meas_df[\"APD_mid_mm\"].notna()].tolist()\nif idx:\n # we need the qc payload from a fresh measurement call (to carry arrays)\n r = perlevel_df.iloc[idx[0]].to_dict()\n m = measure_axial_midline(ct_arr2, ct_img2, seg_arr2, r[\"upper\"], r[\"lower\"], r[\"z_mid\"], r[\"x_mid\"])\n if m and m[\"qc\"]:\n qc_dual_full_and_roi(ct_arr2, ct_img2, seg_arr2, r, m[\"qc\"])\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T22:26:55.805378Z","iopub.execute_input":"2025-09-24T22:26:55.805792Z","iopub.status.idle":"2025-09-24T22:26:57.303766Z","shell.execute_reply.started":"2025-09-24T22:26:55.805768Z","shell.execute_reply":"2025-09-24T22:26:57.302841Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Cell 6H: Robust axial APD/Torg (label∪HU bone, non-bone run) ===\nimport numpy as np\nfrom scipy import ndimage as ndi\n\n# Tunables (adjust if needed)\nHU_BONE_TH = 240.0 # HU threshold to add to bone union (try 220–300 for soft/bone kernels)\nOPEN_SZ = 7 # opening to remove thin posterior elements from vertebra mask\nKEEP_FRAC = 0.65 # keep anterior 65% of vertebra mask as body-only\nZ_NEIGH = 1 # z ∈ [z_mid ± Z_NEIGH]\nX_BAND = 5 # measure also across ±X_BAND cols around midline (robust)\nAREA_MIN = 120 # px: drop tiny bone comps (spicules)\nRUN_LEN = 3 # consecutive bone px to accept \"true\" lamina (de-streak)\nROBUST_Q = 0.75 # percentile for band (robust) APD\nMED_FILT = (3,1) # median filter on ROI to reduce streaks (H×W)\n\ndef _px_from_mm(mm, sp): return int(round(mm / float(sp)))\n\ndef _body_only_mask(ax_seg2d, vid):\n m = (ax_seg2d == int(vid))\n if not m.any(): return None\n opened = ndi.binary_opening(m, structure=np.ones((OPEN_SZ, OPEN_SZ)))\n base = opened if opened.any() else m\n ys = np.where(base)[0]\n y0, y1 = ys.min(), ys.max()\n cut = int(y0 + KEEP_FRAC*(y1-y0))\n body = base.copy(); body[cut+1:, :] = 0\n return body if body.any() else None\n\ndef _refine_midline_x(ax_seg2d, up_id, lo_id, x_default):\n xs=[]\n for vid in (int(up_id), int(lo_id)):\n m=(ax_seg2d==vid)\n if m.any(): xs.append(int(np.median(np.where(m)[1])))\n return int(np.mean(xs)) if xs else int(x_default)\n\ndef _bone_union(ax_seg2d, ax_ct2d):\n bone = (ax_seg2d>0) | (ax_ct2d>HU_BONE_TH)\n # remove tiny comps\n lab, num = ndi.label(bone.astype(np.uint8))\n if num:\n sizes=np.bincount(lab.ravel())\n bone[np.isin(lab, np.where(sizes0:\n apd_mid = gap_px*sy\n vb_ap = max(0.0, (y_post_r - y_ant_r)*sy)\n if vb_ap>0:\n apd_mid_list.append(apd_mid); torg_mid_list.append(apd_mid/vb_ap)\n\n # BAND (robust)\n xb0=max(0, x_mid_r-X_BAND); xb1=min(ct_roi.shape[1]-1, x_mid_r+X_BAND)\n apds=[]; torgs=[]\n for x in range(xb0, xb1+1):\n gpx=_first_nonbone_run_len(nonbone, y_post_r, x)\n if gpx<=0: continue\n apd=gpx*sy; vb_ap=max(0.0, (y_post_r-y_ant_r)*sy)\n if vb_ap>0: apds.append(apd); torgs.append(apd/vb_ap)\n if apds:\n apd_band_list.append(float(np.quantile(apds, ROBUST_Q)))\n torg_band_list.append(float(np.quantile(torgs, ROBUST_Q)))\n\n qc_payload=dict(z=z, y0=y0, y1=y1, x0=x0, x1=x1, x_mid_r=x_mid_r,\n y_ant_r=y_ant_r, y_post_r=y_post_r, ct_roi=ct_roi,\n seg_roi=seg_roi, nonbone=nonbone)\n\n if not (apd_mid_list or apd_band_list): return None\n return {\n \"APD_mid_mm\": np.median(apd_mid_list) if apd_mid_list else np.nan,\n \"APD_band_mm\": np.median(apd_band_list) if apd_band_list else np.nan,\n \"Torg_mid\": np.median(torg_mid_list) if torg_mid_list else np.nan,\n \"Torg_band\": np.median(torg_band_list)if torg_band_list else np.nan,\n \"z_mid_used\": z_mid, \"x_mid_used\": x_mid, \"qc\": qc_payload\n }\n\ndef sev_from_apd(a): \n if np.isnan(a): return -1\n return 2 if a<10 else (1 if a<13 else 0)\n\n# ---- recompute on current study using robust method\nrows=[]\nfor r in perlevel_df.to_dict(\"records\"):\n m = measure_axial_robust(ct_arr2, ct_img2, seg_arr2, r[\"upper\"], r[\"lower\"], r[\"z_mid\"], r[\"x_mid\"])\n if m is None:\n rows.append({**r, \"APD_mid_mm\":np.nan, \"APD_band_mm\":np.nan, \"Torg_mid\":np.nan, \"Torg_band\":np.nan,\n \"severity_mid\":-1, \"severity_band\":-1})\n else:\n rows.append({**r, **{k:v for k,v in m.items() if k!=\"qc\"},\n \"severity_mid\": sev_from_apd(m[\"APD_mid_mm\"]),\n \"severity_band\": sev_from_apd(m[\"APD_band_mm\"])})\nmeas_df = pd.DataFrame(rows)\nprint(\"[MEAS-robust] APD (midline & band) and severities:\")\ndisplay(meas_df)\n\n# ---- dual-view QC (full + ROI) reusing 6G but with robust non-bone overlay\ndef qc_dual_full_and_roi_robust(ct_arr, ct_img, seg_arr, level_row, qc, window_full=(400,1800), window_roi=(500,2200)):\n sx, sy, sz = ct_img.GetSpacing()\n z=qc[\"z\"]; y0,y1,x0,x1 = qc[\"y0\"],qc[\"y1\"],qc[\"x0\"],qc[\"x1\"]\n x_mid_r=qc[\"x_mid_r\"]; y_ant_r=qc[\"y_ant_r\"]; y_post_r=qc[\"y_post_r\"]\n ct_roi=qc[\"ct_roi\"]; seg_roi=qc[\"seg_roi\"]; nonbone=qc[\"nonbone\"]\n\n # full view\n full=window_image(ct_arr[z], *window_full); seg_full=seg_arr[z]\n fig,axs=plt.subplots(1,2,figsize=(11,5.5))\n extent_full=[0, full.shape[1]*sx, full.shape[0]*sy, 0]\n axs[0].imshow(overlay_segmentation(full, seg_full), extent=extent_full); axs[0].axis('off'); axs[0].set_aspect('equal')\n axs[0].axvline((x0+x_mid_r)*sx, color='white', lw=1.2, ls='--')\n axs[0].set_title(f\"Full axial (z={z})\")\n\n # ROI\n roi=window_image(ct_roi, *window_roi)\n extent_roi=[0, roi.shape[1]*sx, roi.shape[0]*sy, 0]\n axs[1].imshow(overlay_segmentation(roi, seg_roi), extent=extent_roi); axs[1].axis('off'); axs[1].set_aspect('equal')\n # overlay non-bone (canal) as faint green\n canal = np.ma.masked_where(nonbone==0, nonbone)\n axs[1].imshow(canal, extent=extent_roi, cmap='Greens', alpha=0.25, interpolation='nearest')\n x_mm=x_mid_r*sx\n axs[1].axvline(x_mm, color='white', lw=1.2, ls='--')\n axs[1].hlines([y_ant_r*sy, y_post_r*sy], x_mm-8, x_mm+8, colors=['cyan','yellow'], lw=2)\n # compute midline gap for display\n gap_px=_first_nonbone_run_len(nonbone, y_post_r, x_mid_r)\n if gap_px>0:\n apd=gap_px*sy\n axs[1].hlines([(y_post_r+gap_px)*sy], x_mm-8, x_mm+8, colors='red', lw=2)\n axs[1].text(x_mm+6, (y_post_r*sy + (y_post_r+gap_px)*sy)/2,\n f\"APD(mid) ≈ {apd:.1f} mm\", color='white',\n bbox=dict(fc='black', alpha=0.55), fontsize=10)\n axs[1].set_title(\"Coned-in ROI (robust canal)\")\n\n plt.suptitle(f\"QC — {level_row['level']} (midline & ROI)\", y=0.98, fontsize=12)\n plt.tight_layout(); plt.show()\n\n# Show QC for first measurable level\nidx = meas_df.index[meas_df[\"APD_mid_mm\"].notna()].tolist()\nif idx:\n r = perlevel_df.iloc[idx[0]].to_dict()\n m = measure_axial_robust(ct_arr2, ct_img2, seg_arr2, r[\"upper\"], r[\"lower\"], r[\"z_mid\"], r[\"x_mid\"])\n if m and m[\"qc\"]:\n qc_dual_full_and_roi_robust(ct_arr2, ct_img2, seg_arr2, r, m[\"qc\"])\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T22:30:13.88006Z","iopub.execute_input":"2025-09-24T22:30:13.880921Z","iopub.status.idle":"2025-09-24T22:30:15.534201Z","shell.execute_reply.started":"2025-09-24T22:30:13.880891Z","shell.execute_reply":"2025-09-24T22:30:15.533299Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Cell 7H: Batch runner for all eligible studies with cervical segs ===\nfrom tqdm.auto import tqdm\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\n# Build \"eligible\" list once (has seg and at least one C1..C7 voxel)\neligible = []\nfor p in study_dirs:\n if p.name not in seg_map: \n continue\n img, arr = load_seg_nifti(seg_map[p.name])\n if np.any(np.isin(arr, np.arange(1,8))):\n eligible.append(p)\nprint(f\"[BATCH] Eligible cervical studies: {len(eligible)}\")\n\n# Process N studies (set N=None for all)\nN = None\nrows_all = []\nfor p in tqdm(eligible[:N], desc=\"Measuring APD\"):\n # load best CT series\n ct_img, ct_arr, files, uid = load_best_series_any(p)\n if ct_img is None: \n continue\n seg_img, seg_arr = load_seg_nifti(seg_map[p.name])\n # resample seg→CT if needed\n same_geom = (ct_img.GetSize()==seg_img.GetSize() and\n np.allclose(ct_img.GetSpacing(), seg_img.GetSpacing(), atol=1e-3) and\n tuple(ct_img.GetDirection())==tuple(seg_img.GetDirection()))\n if not same_geom:\n seg_arr = sitk.GetArrayFromImage(resample_like(seg_img, ct_img, True))\n\n # per-level planes (C2/3..C6/7)\n recs=[]\n for lid in range(2,8):\n upper, lower = lid, min(lid+1,7)\n # disc mid-z\n up = (seg_arr==upper); lo=(seg_arr==lower)\n z_mid = int(np.median(np.argwhere(up)[:,0])) if up.any() else (int(np.median(np.argwhere(lo)[:,0])) if lo.any() else seg_arr.shape[0]//2)\n # initial x at upper body center; refined inside measure_axial_robust\n x_mid = int(np.median(np.where(up)[1])) if up.any() else seg_arr.shape[2]//2\n recs.append({\"study_id\": p.name, \"level\": f\"C{upper}/C{lower}\", \"upper\":upper, \"lower\":lower,\n \"z_mid\": z_mid, \"x_mid\": x_mid})\n perlevel_df_batch = pd.DataFrame(recs)\n\n # robust measurements\n for r in perlevel_df_batch.to_dict(\"records\"):\n m = measure_axial_robust(ct_arr, ct_img, seg_arr, r[\"upper\"], r[\"lower\"], r[\"z_mid\"], r[\"x_mid\"])\n if m is None:\n rows_all.append({**r, \"APD_mid_mm\":np.nan, \"APD_band_mm\":np.nan, \"Torg_mid\":np.nan, \"Torg_band\":np.nan,\n \"severity_mid\":-1, \"severity_band\":-1})\n else:\n rows_all.append({**r, **{k:v for k,v in m.items() if k!=\"qc\"},\n \"severity_mid\": sev_from_apd(m[\"APD_mid_mm\"]),\n \"severity_band\": sev_from_apd(m[\"APD_band_mm\"])})\n# Save cohort CSV\ncohort_df = pd.DataFrame(rows_all)\ncsv_out = \"/kaggle/working/ccs_cohort_apd.csv\"\ncohort_df.to_csv(csv_out, index=False)\nprint(f\"[SAVE] Cohort measurements -> {csv_out}\\nshape={cohort_df.shape}\")\n\n# Quick plots (overall)\nok_mid = cohort_df[\"APD_mid_mm\"].dropna()\nok_band = cohort_df[\"APD_band_mm\"].dropna()\nplt.figure(figsize=(7,3)); plt.hist(ok_band, bins=30, alpha=0.8)\nplt.axvline(10, color=\"r\", ls=\"--\", label=\"Severe < 10 mm\"); plt.axvline(13, color=\"orange\", ls=\"--\", label=\"Moderate 10–13 mm\")\nplt.xlabel(\"AP canal diameter (mm), band-robust\"); plt.ylabel(\"Count\"); plt.legend(); plt.title(\"Cervical APD (band-robust) — cohort\")\nplt.tight_layout(); plt.show()\n\n# Per-level medians\nmed = cohort_df.groupby(\"level\")[\"APD_band_mm\"].median().reset_index()\nprint(\"[SUMMARY] Median APD_band_mm by level:\"); display(med)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T22:30:45.960625Z","iopub.execute_input":"2025-09-24T22:30:45.961225Z","iopub.status.idle":"2025-09-24T22:44:16.541977Z","shell.execute_reply.started":"2025-09-24T22:30:45.961198Z","shell.execute_reply":"2025-09-24T22:44:16.539531Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Final: fast batch + report (no training) ===\nimport numpy as np, pandas as pd, matplotlib.pyplot as plt\nfrom tqdm.auto import tqdm\n\n# ---------- knobs ----------\nFAST_N = 10 # how many eligible studies to process quickly\nLEVEL_WHITELIST = {\"C4/C5\",\"C5/C6\",\"C6/C7\"} # focus cohort metrics on these\nQC_EXAMPLES = 3 # number of dual-view QC panels to render\nRANDOM_SEED = 42\nnp.random.seed(RANDOM_SEED)\n\n# ---------- safety checks ----------\nassert 'study_dirs' in globals() and 'seg_map' in globals(), \"Run the EDA cells first.\"\nassert 'load_seg_nifti' in globals() and 'load_best_series_any' in globals(), \"Run the loader cells.\"\nassert 'resample_like' in globals(), \"Missing resample_like.\"\nassert 'measure_axial_robust' in globals(), \"Run the robust measurement cell (6H).\"\n\n# ---------- discover eligible (has segmentation and some C1..C7 voxels) ----------\neligible = []\nfor p in study_dirs:\n if p.name not in seg_map:\n continue\n _, seg_arr0 = load_seg_nifti(seg_map[p.name])\n if np.any(np.isin(seg_arr0, np.arange(1,8))):\n eligible.append(p)\nprint(f\"[BATCH] Eligible cervical studies with segs: {len(eligible)}\")\n\n# ---------- process a fast subset ----------\nrows_all = []\nuse_list = eligible[:FAST_N] if FAST_N is not None else eligible\nprint(f\"[RUN] Measuring {len(use_list)} studies...\")\n\nfor p in tqdm(use_list, desc=\"APD measure\"):\n # 1) CT series\n ct_img, ct_arr, files, uid = load_best_series_any(p)\n if ct_img is None:\n continue\n # 2) segmentation (NIfTI)\n seg_img, seg_arr = load_seg_nifti(seg_map[p.name])\n # 3) resample seg -> CT geometry if needed\n same_geom = (ct_img.GetSize()==seg_img.GetSize() and\n np.allclose(ct_img.GetSpacing(), seg_img.GetSpacing(), atol=1e-3) and\n tuple(ct_img.GetDirection())==tuple(seg_img.GetDirection()))\n if not same_geom:\n seg_arr = sitk.GetArrayFromImage(resample_like(seg_img, ct_img, True))\n\n # 4) define disc planes C2/3..C6/7 (simple, robust)\n perlevel = []\n for lid in range(2,8): # C2..C7 => planes C2/3..C7/7 (C7/7 ignored later)\n upper, lower = lid, min(lid+1,7)\n up = (seg_arr==upper); lo=(seg_arr==lower)\n if up.any():\n z_mid = int(np.median(np.argwhere(up)[:,0]))\n x_mid = int(np.median(np.where(up)[1]))\n elif lo.any():\n z_mid = int(np.median(np.argwhere(lo)[:,0]))\n x_mid = int(np.median(np.where(lo)[1]))\n else:\n z_mid = seg_arr.shape[0]//2\n x_mid = seg_arr.shape[2]//2\n perlevel.append({\"study_id\": p.name, \"level\": f\"C{upper}/C{lower}\",\n \"upper\":upper, \"lower\":lower, \"z_mid\":z_mid, \"x_mid\":x_mid})\n perlevel_df_batch = pd.DataFrame(perlevel)\n\n # 5) robust measurements\n for r in perlevel_df_batch.to_dict(\"records\"):\n m = measure_axial_robust(ct_arr, ct_img, seg_arr, r[\"upper\"], r[\"lower\"], r[\"z_mid\"], r[\"x_mid\"])\n if m is None:\n rows_all.append({**r, \"APD_mid_mm\":np.nan, \"APD_band_mm\":np.nan, \"Torg_mid\":np.nan, \"Torg_band\":np.nan})\n else:\n rows_all.append({**r, **{k:v for k,v in m.items() if k!=\"qc\"}})\n\n# ---------- tidy & final severity ----------\ndef sev_from_apd(a):\n if np.isnan(a): return -1\n return 2 if a < 10.0 else (1 if a < 13.0 else 0)\n\ncohort_df = pd.DataFrame(rows_all)\n# choose band-robust when available, else midline\ncohort_df[\"APD_final_mm\"] = np.where(cohort_df[\"APD_band_mm\"].notna(),\n cohort_df[\"APD_band_mm\"], cohort_df[\"APD_mid_mm\"])\ncohort_df[\"severity_final\"] = cohort_df[\"APD_final_mm\"].apply(sev_from_apd)\n\ncsv_out = \"/kaggle/working/ccs_cohort_apd_final.csv\"\ncohort_df.to_csv(csv_out, index=False)\nprint(f\"[SAVE] Cohort CSV -> {csv_out} shape={cohort_df.shape}\")\n\n# ---------- headline summary (focus on C4–C7) ----------\ncohort_eval = cohort_df[cohort_df[\"level\"].isin(LEVEL_WHITELIST)].copy()\nn_levels = cohort_eval[\"APD_final_mm\"].notna().sum()\nn_studies = cohort_eval[\"study_id\"].nunique()\ntriage_flag = (cohort_eval\n .assign(flag=lambda d: d[\"severity_final\"].ge(1))\n .groupby(\"study_id\")[\"flag\"].any()\n .mean())\nprint(f\"[SUMMARY] Levels with measurements (C4–C7): {n_levels} across {n_studies} studies\")\nprint(f\"[SUMMARY] % patients flagged (moderate+severe at any C4–C7): {100*triage_flag:.1f}%\")\n\nmedians = (cohort_eval\n .groupby(\"level\")[\"APD_final_mm\"]\n .median()\n .reindex(sorted(LEVEL_WHITELIST)))\nprint(\"[SUMMARY] Per-level median APD_final_mm (mm):\")\nprint(medians)\n\n# ---------- quick plots ----------\nok = cohort_eval[\"APD_final_mm\"].dropna()\nplt.figure(figsize=(7,3.3))\nplt.hist(ok, bins=28, alpha=0.85)\nplt.axvline(10, color=\"r\", ls=\"--\", label=\"Severe < 10 mm\")\nplt.axvline(13, color=\"orange\", ls=\"--\", label=\"Moderate 10–13 mm\")\nplt.xlabel(\"AP canal diameter (mm) — APD_final\"); plt.ylabel(\"Count\"); plt.legend()\nplt.title(\"Cervical APD (C4–C7) — cohort\")\nplt.tight_layout(); plt.show()\n\nplt.figure(figsize=(6,3.3))\nlvls = sorted(list(LEVEL_WHITELIST))\ndata = [cohort_eval.loc[cohort_eval[\"level\"]==L, \"APD_final_mm\"].dropna().values for L in lvls]\nplt.boxplot(data, labels=lvls, showmeans=True)\nplt.ylabel(\"AP canal diameter (mm)\"); plt.title(\"Per-level APD (C4–C7)\")\nplt.tight_layout(); plt.show()\n\n# ---------- (optional) QC dual-view panels for 3 random measured levels ----------\nif QC_EXAMPLES > 0:\n # helper (re-uses the robust function to get QC payload)\n def qc_dual_full_and_roi_robust(ct_arr, ct_img, seg_arr, level_row, qc,\n window_full=(400,1800), window_roi=(500,2200)):\n sx, sy, sz = ct_img.GetSpacing()\n z=qc[\"z\"]; y0,y1,x0,x1 = qc[\"y0\"],qc[\"y1\"],qc[\"x0\"],qc[\"x1\"]\n x_mid_r=qc[\"x_mid_r\"]; y_ant_r=qc[\"y_ant_r\"]; y_post_r=qc[\"y_post_r\"]\n ct_roi=qc[\"ct_roi\"]; seg_roi=qc[\"seg_roi\"]; nonbone=qc[\"nonbone\"]\n\n # full-FOV\n full=window_image(ct_arr[z], *window_full); seg_full=seg_arr[z]\n fig,axs=plt.subplots(1,2,figsize=(11,5.5))\n extent_full=[0, full.shape[1]*sx, full.shape[0]*sy, 0]\n axs[0].imshow(overlay_segmentation(full, seg_full), extent=extent_full); axs[0].axis('off'); axs[0].set_aspect('equal')\n axs[0].axvline((x0+x_mid_r)*sx, color='white', lw=1.2, ls='--'); axs[0].set_title(f\"Full axial (z={z})\")\n\n # ROI\n roi=window_image(ct_roi, *window_roi); extent_roi=[0, roi.shape[1]*sx, roi.shape[0]*sy, 0]\n axs[1].imshow(overlay_segmentation(roi, seg_roi), extent=extent_roi); axs[1].axis('off'); axs[1].set_aspect('equal')\n canal = np.ma.masked_where(nonbone==0, nonbone)\n axs[1].imshow(canal, extent=extent_roi, cmap='Greens', alpha=0.25, interpolation='nearest')\n x_mm=x_mid_r*sx\n axs[1].axvline(x_mm, color='white', lw=1.2, ls='--')\n axs[1].hlines([y_ant_r*sy, y_post_r*sy], x_mm-8, x_mm+8, colors=['cyan','yellow'], lw=2)\n # midline gap for display\n gap_px = 0\n H = nonbone.shape[0]\n for yy in range(int(y_post_r)+1, H):\n if nonbone[yy, int(x_mid_r)]: gap_px += 1\n else: break\n if gap_px>0:\n apd=gap_px*sy\n axs[1].hlines([(y_post_r+gap_px)*sy], x_mm-8, x_mm+8, colors='red', lw=2)\n axs[1].text(x_mm+6, (y_post_r*sy + (y_post_r+gap_px)*sy)/2,\n f\"APD(mid) ≈ {apd:.1f} mm\", color='white',\n bbox=dict(fc='black', alpha=0.55), fontsize=10)\n axs[1].set_title(\"Coned-in ROI (robust canal)\")\n plt.suptitle(f\"QC — {level_row['level']} (midline & ROI)\", y=0.98, fontsize=12)\n plt.tight_layout(); plt.show()\n\n # sample levels with measurements\n candidates = cohort_eval.dropna(subset=[\"APD_final_mm\"]).sample(min(QC_EXAMPLES, \n cohort_eval[\"APD_final_mm\"].notna().sum()), random_state=RANDOM_SEED)\n print(f\"[QC] Rendering {len(candidates)} dual-view panels...\")\n for _, rr in candidates.iterrows():\n sid = rr[\"study_id\"]; Ltxt = rr[\"level\"]; upper=int(Ltxt.split('/')[0][1]); lower=int(Ltxt.split('/')[1][1])\n # reload study (simple and reliable for small QC count)\n p = [pp for pp in eligible if pp.name==sid][0]\n ct_img, ct_arr, _, _ = load_best_series_any(p)\n seg_img, seg_arr = load_seg_nifti(seg_map[p.name])\n same_geom = (ct_img.GetSize()==seg_img.GetSize() and\n np.allclose(ct_img.GetSpacing(), seg_img.GetSpacing(), atol=1e-3) and\n tuple(ct_img.GetDirection())==tuple(seg_img.GetDirection()))\n if not same_geom:\n seg_arr = sitk.GetArrayFromImage(resample_like(seg_img, ct_img, True))\n r0 = {\"upper\":upper, \"lower\":lower, \"z_mid\":int(rr[\"z_mid_used\"]), \"x_mid\":int(rr[\"x_mid_used\"])}\n m = measure_axial_robust(ct_arr, ct_img, seg_arr, **r0)\n if m and m[\"qc\"]:\n r_level = {\"level\": Ltxt}\n qc_dual_full_and_roi_robust(ct_arr, ct_img, seg_arr, r_level, m[\"qc\"])\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Notes / citations: For clinically interpretable thresholds, we use AP canal diameter at the disc plane with severe < 10 mm and moderate 10–13 mm, plus Torg–Pavlov ratio < 0.8 as supportive criteria—standard in cervical stenosis literature and appropriate for CT-based triage in MRI-limited settings (Pavlov & Torg, 1987). The RSNA-2022 Cervical Spine CT dataset supplies vertebra segmentations (C1–C7/T1–T12) but no stenosis ground truth; this pipeline therefore provides deterministic morphometrics and cohort summaries that can serve as silver labels and as an operational triage pilot.","metadata":{}},{"cell_type":"code","source":"# === Cell: MRI → CT transfer weights resolver (no training) ===\nimport os, glob, torch, timm\n\n# ---- (A) optional: set your known MRI checkpoint path here (if you added a dataset with weights)\nPREFER_USER_CKPT = \"/kaggle/input/lumbar-model-weights/_scs_classify_5ch_axsagt2-lstm-mil_auxloss_auxdepth_convnext-s_for_exp0.ckpt\"\n# Leave as None to auto-search\nif not os.path.exists(PREFER_USER_CKPT):\n PREFER_USER_CKPT = None\n\n# ---- (B) auto-search common locations if not specified\nSEARCH_GLOBS = [\n \"/kaggle/input/*mri*weights*/*.pt\",\n \"/kaggle/input/*mri*weights*/*.pth\",\n \"/kaggle/input/*lumbar*weights*/*.pt\",\n \"/kaggle/input/*lumbar*weights*/*.pth\",\n \"/kaggle/input/*weights*/*.ckpt\",\n \"/kaggle/input/*mri*/*.ckpt\",\n]\ndef find_mri_ckpt():\n if PREFER_USER_CKPT and os.path.exists(PREFER_USER_CKPT):\n return PREFER_USER_CKPT\n for pat in SEARCH_GLOBS:\n hits = sorted(glob.glob(pat))\n if hits:\n return hits[0]\n return None\n\n# ---- (C) robust loader that adapts first conv to 1-channel and strips Lightning prefixes\ndef load_mri_transfer_into_backbone(backbone, ckpt_path, in_chans_expected=1):\n if ckpt_path is None or not os.path.exists(ckpt_path):\n print(\"[xfer] No MRI checkpoint found; using ImageNet weights already in timm.\")\n return 0\n sd = torch.load(ckpt_path, map_location=\"cpu\")\n if \"state_dict\" in sd:\n sd = sd[\"state_dict\"]\n # Keep only backbone.* if present; strip prefix\n cleaned = {}\n for k,v in sd.items():\n if \"backbone\" in k:\n k2 = k.split(\"backbone.\",1)[-1]\n cleaned[k2] = v\n elif k in backbone.state_dict():\n cleaned[k] = v\n\n # First conv channel adaptation if needed (e.g., 5ch → 1ch)\n for k,v in list(cleaned.items()):\n if \"weight\" in k and v.ndim==4 and v.shape[1] != in_chans_expected:\n with torch.no_grad():\n if v.shape[1] > in_chans_expected:\n v = v.mean(dim=1, keepdim=True).repeat(1, in_chans_expected, 1, 1)\n else:\n v = v.repeat(1, in_chans_expected//v.shape[1] + 1, 1, 1)[:, :in_chans_expected]\n cleaned[k] = v\n\n missing, unexpected = backbone.load_state_dict(cleaned, strict=False)\n loaded = len(cleaned) - len(missing)\n print(f\"[xfer] Loaded {loaded} tensors from MRI ckpt: {os.path.basename(ckpt_path)} \"\n f\"(missing={len(missing)}, unexpected={len(unexpected)})\")\n return loaded\n\n# ---- (D) build backbone and apply transfer (or ImageNet)\nMODEL_NAME = \"convnext_tiny.in12k_ft_in1k\" # good small backbone\nIN_CHANS = 1 # we use single-channel axial slices\nprint(f\"[model] Building {MODEL_NAME} (ImageNet weights) with in_chans={IN_CHANS}\")\nbackbone = timm.create_model(MODEL_NAME, pretrained=True, in_chans=IN_CHANS, num_classes=0)\n\nMRI_CKPT = find_mri_ckpt()\nif MRI_CKPT:\n _ = load_mri_transfer_into_backbone(backbone, MRI_CKPT, in_chans_expected=IN_CHANS)\nelse:\n print(\"[xfer] MRI weights not found; proceeding with ImageNet initialization.\")\n\n# ---- (E) tiny sanity: forward a dummy slice to confirm the feature size\nwith torch.no_grad():\n dummy = torch.randn(1, IN_CHANS, 256, 256)\n feat = backbone(dummy)\n if feat.ndim == 4:\n feat = torch.nn.AdaptiveAvgPool2d(1)(feat).flatten(1)\n nf = feat.shape[1]\nprint(f\"[model] Backbone feature dimension: {nf}\")\n# Save the backbone state_dict so downstream cells can reuse it quickly if needed\ntorch.save(backbone.state_dict(), \"/kaggle/working/backbone_ct_init.pth\")\nprint(\"[save] /kaggle/working/backbone_ct_init.pth\")\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# After you have `meas_df` from 6H/6G\nimport numpy as np\n\nmeas_df[\"APD_final_mm\"] = np.where(\n meas_df[\"APD_band_mm\"].notna(), meas_df[\"APD_band_mm\"], meas_df[\"APD_mid_mm\"]\n)\n\ndef sev_from_apd(a):\n if np.isnan(a): return -1\n return 2 if a < 10.0 else (1 if a < 13.0 else 0) # severe <10, moderate 10–13\nmeas_df[\"severity_final\"] = meas_df[\"APD_final_mm\"].apply(sev_from_apd)\n\n# OPTIONAL clamp: if APD_mid << APD_band by >6 mm (likely artifact), trust the band\nmeas_df[\"APD_final_mm\"] = np.where(\n (meas_df[\"APD_mid_mm\"].notna()) & (meas_df[\"APD_band_mm\"].notna()) &\n ((meas_df[\"APD_band_mm\"] - meas_df[\"APD_mid_mm\"]) > 6.0),\n meas_df[\"APD_band_mm\"], meas_df[\"APD_final_mm\"]\n)\nmeas_df[\"severity_final\"] = meas_df[\"APD_final_mm\"].apply(sev_from_apd)\ndisplay(meas_df[[\"level\",\"APD_mid_mm\",\"APD_band_mm\",\"APD_final_mm\",\"severity_mid\",\"severity_band\",\"severity_final\"]])\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Cell 7A: QC overlay for a chosen level (draw APD) ===\ndef qc_draw_level(ct_arr, ct_img, seg_arr, level_row, flip_ud=True):\n lvl = level_row[\"level\"]; z_mid, x_mid = int(level_row[\"z_mid\"]), int(level_row[\"x_mid\"])\n sx, sy, sz = ct_img.GetSpacing()\n # build a sagittal band around z_mid for nicer context\n band = ct_arr[max(0,z_mid-2):min(ct_arr.shape[0], z_mid+3), :, x_mid].astype(np.float32)\n # collapse band to median for display\n sag_disp = np.median(band, axis=0) # (Y,)\n # expand to 2D image for overlay (H, W) as (Zband, Y)\n sag2d = band # (Hk, Y)\n # window & overlay vertebra labels projected to this slice column\n lbl_band = seg_arr[max(0,z_mid-2):min(ct_arr.shape[0], z_mid+3), :, x_mid]\n lbl_proj = (lbl_band.max(axis=0)).astype(np.int32) # (Y,)\n\n # Construct 2D for overlay function: make HxW image (use H=band size)\n v = np.tile(sag_disp[None,:], (lbl_band.shape[0],1))\n ov = overlay_segmentation(window_image(v, 400, 1800), np.tile(lbl_proj[None,:], (lbl_band.shape[0],1)))\n\n # Compute edges for a single z within band (center slice)\n # inside qc_draw_level(...)\n res_vb = _posterior_body_from_mask(lbl_band, level_row[\"upper\"])\n y_lam_ant = _lamina_anterior_y(band, res_vb[1]) if res_vb is not None else None\n\n fig, ax = plt.subplots(figsize=(6, 6 * (lbl_band.shape[0]*sz)/(lbl_proj.shape[0]*sy)))\n extent = [0, lbl_proj.shape[0]*sy, lbl_band.shape[0]*sz, 0]\n ax.imshow(ov, extent=extent)\n if res is not None:\n y_vb_ant, y_vb_post, y_lam_ant = res\n # draw vertical lines (y axis in image coordinate)\n yv = [y_vb_ant*sy, y_vb_post*sy, y_lam_ant*sy]\n for i,(c,lbl) in enumerate(zip([\"cyan\",\"yellow\",\"red\"], [\"VB anterior\",\"VB posterior\",\"Lamina anterior\"])):\n ax.axvline(yv[i], color=c, lw=2, label=lbl)\n apd_mm = max(0.0, yv[2]-yv[1])\n ax.text(5, 5, f\"APD ≈ {apd_mm:.1f} mm\", color=\"white\", fontsize=11, ha=\"left\", va=\"top\",\n bbox=dict(facecolor=\"black\", alpha=0.5))\n ax.legend(loc=\"lower right\", frameon=False)\n ax.set_title(f\"QC — {lvl} (sagittal band at x={x_mid})\")\n ax.set_xlabel(\"mm (anterior → posterior)\"); ax.set_ylabel(\"mm (superior → inferior)\")\n ax.axis('off'); plt.tight_layout(); plt.show()\n\n# Draw QC for the first level that has a measurement\nok_rows = meas_df[meas_df[\"severity\"] >= 0]\nif len(ok_rows):\n qc_draw_level(ct_arr2, ct_img2, seg_arr2, ok_rows.iloc[0])\nelse:\n print(\"No measurable levels in this study.\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T22:24:53.477875Z","iopub.execute_input":"2025-09-24T22:24:53.47822Z","iopub.status.idle":"2025-09-24T22:24:53.628754Z","shell.execute_reply.started":"2025-09-24T22:24:53.4782Z","shell.execute_reply":"2025-09-24T22:24:53.627483Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Cell 8A: Export + summary ===\nOUT_DF = pd.DataFrame(meas_rows)\nOUT_DF[\"study_id\"] = study_dir2.name\ncols = [\"study_id\",\"level\",\"APD_mm\",\"Torg\",\"severity\",\"z_mid\",\"x_mid\",\"upper\",\"lower\"]\nOUT_DF = OUT_DF[cols]\ncsv_path = f\"/kaggle/working/ccs_measurements_{study_dir2.name}.csv\"\nOUT_DF.to_csv(csv_path, index=False)\nprint(f\"[SAVE] Wrote per-level measurements: {csv_path}\")\n\n# Summary bar (severity per level for this study)\nsev_map = {0:\"None/Mild\",1:\"Moderate\",2:\"Severe\",-1:\"N/A\"}\nOUT_DF[\"severity_txt\"] = OUT_DF[\"severity\"].map(sev_map)\nax = OUT_DF.plot(x=\"level\", y=\"APD_mm\", kind=\"bar\", legend=False, figsize=(7,3.2))\nax.set_ylim(0, 25); ax.set_ylabel(\"AP canal diameter (mm)\"); ax.set_xlabel(\"Level\")\nfor p, s in zip(ax.patches, OUT_DF[\"severity_txt\"]):\n ax.annotate(s, (p.get_x()+p.get_width()/2, p.get_height()+0.3), ha=\"center\", fontsize=9)\nplt.tight_layout(); plt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T22:24:56.181178Z","iopub.execute_input":"2025-09-24T22:24:56.182045Z","iopub.status.idle":"2025-09-24T22:24:56.370929Z","shell.execute_reply.started":"2025-09-24T22:24:56.182017Z","shell.execute_reply":"2025-09-24T22:24:56.370127Z"}},"outputs":[],"execution_count":null}]}