From 476f395c5a6bc6618e57f705ca36ba5524a44d69 Mon Sep 17 00:00:00 2001 From: Petra Hedesiu Date: Wed, 27 May 2026 08:47:01 -0400 Subject: [PATCH 1/5] Add post-hoc AEA classification with bone-contact roof detection Classifies 4 binary clinical labels (roof contact L/R, ethmoid filled L/R) from segmentation masks using logistic regression. Roof contact detection uses CT bone HU sampling above the mask superior surface (bone_contact_ratio, mean_above_hu) instead of z-position alone, lifting AUC from ~0.42 to ~0.75. Co-Authored-By: Claude Sonnet 4.6 --- try2/experiments/exp_classification.py | 335 +++++++++++++++++++++++++ try2/shared/classification.py | 264 +++++++++++++++++++ 2 files changed, 599 insertions(+) create mode 100644 try2/experiments/exp_classification.py create mode 100644 try2/shared/classification.py diff --git a/try2/experiments/exp_classification.py b/try2/experiments/exp_classification.py new file mode 100644 index 0000000..84ad008 --- /dev/null +++ b/try2/experiments/exp_classification.py @@ -0,0 +1,335 @@ +""" +Train and evaluate AEA classification from segmentation masks. + +Usage: + # Train + LOO-CV on CROP1, then optionally validate on CROP EXT: + python exp_classification.py [--data-dir DATA_DIR] [--xlsx PATH] [--output-dir DIR] + [--ext-dir EXT_DIR] + +CROP1 folder format : "016. VA016 VINTELER ANA-MARIA" → code extracted as "VA016" +CROP EXT structure : + EXT_DIR/ + 01/ + EXT01/ ← DICOM directory + EXT01.nrrd ← segmentation mask + 02/ + EXT02/ + EXT02.nrrd + ... + +EXT patients are matched to the Excel by their EXTxx code. If the Excel has no +EXT codes, predictions are still written to output_dir/ext_predictions.csv. +""" + +import sys +import os +import argparse +import re +import csv + +sys.path.insert(0, os.path.join(os.path.dirname(__file__), '..', '..')) +sys.path.insert(0, os.path.join(os.path.dirname(__file__), '..')) + +import numpy as np +from tqdm import tqdm +from sklearn.model_selection import LeaveOneOut +from sklearn.metrics import ( + accuracy_score, roc_auc_score, classification_report, balanced_accuracy_score +) + +from data_utils import load_patient_data, discover_patients, get_labeled_slice_indices +from shared.classification import ( + load_classification_labels, + build_feature_matrix, + train_classifiers, + save_classifiers, + predict_classifications, + TASK_TO_FEATURES, + _make_clf, +) +from shared.config import ExperimentConfig + + +XLSX_PATH = os.path.join( + os.path.dirname(__file__), '..', '..', 'AEA_classification_sinuzite.xlsx' +) + + +# --------------------------------------------------------------------------- +# Patient code extraction +# --------------------------------------------------------------------------- + +def _extract_patient_code(dicom_dir: str) -> str: + """Folder format: "016. VA016 VINTELER ANA-MARIA" — search anywhere in name.""" + basename = os.path.basename(os.path.normpath(dicom_dir)) + m = re.search(r'([A-Z]{2}\d{3})', basename.upper()) + if m: + return m.group(1) + return basename.upper().strip() + + +# --------------------------------------------------------------------------- +# CROP1 data loading +# --------------------------------------------------------------------------- + +def load_data(data_dir: str): + patients_meta = discover_patients(data_dir) + volumes, masks, codes = [], [], [] + for p in tqdm(patients_meta, desc="Loading CROP1"): + try: + vol, seg, meta = load_patient_data( + p['dicom_dir'], p['nrrd_path'], verbose=False + ) + except Exception as e: + print(f" [skip] {p.get('dicom_dir','?')}: {e}") + continue + if not meta.get('alignment_success', False): + continue + if len(get_labeled_slice_indices(seg)) < 2: + continue + volumes.append(vol) + masks.append(seg) + codes.append(_extract_patient_code(p['dicom_dir'])) + return volumes, masks, codes + + +# --------------------------------------------------------------------------- +# CROP EXT data loading +# --------------------------------------------------------------------------- + +def discover_ext_patients(ext_dir: str) -> list: + """ + Walk CROP EXT structure: + ext_dir/01/EXT01/ (DICOM) + ext_dir/01/EXT01.nrrd + ext_dir/02/EXT02/ + ext_dir/02/EXT02.nrrd + ... + Returns list of dicts: {code, dicom_dir, nrrd_path} + """ + patients = [] + for index_folder in sorted(os.listdir(ext_dir)): + index_path = os.path.join(ext_dir, index_folder) + if not os.path.isdir(index_path): + continue + + # Find the EXTxx subfolder and matching nrrd + ext_code = None + dicom_dir = None + nrrd_path = None + + for entry in os.listdir(index_path): + entry_path = os.path.join(index_path, entry) + if os.path.isdir(entry_path) and re.match(r'^EXT\d+$', entry.upper()): + dicom_dir = entry_path + ext_code = entry.upper() + elif entry.lower().endswith('.nrrd') and re.match(r'^EXT\d+', entry.upper()): + nrrd_path = entry_path + + if dicom_dir and nrrd_path and ext_code: + patients.append({ + 'code': ext_code, + 'dicom_dir': dicom_dir, + 'nrrd_path': nrrd_path, + }) + else: + print(f" [skip EXT] {index_path}: could not find EXTxx dir + nrrd") + + return patients + + +def load_ext_data(ext_dir: str): + patients_meta = discover_ext_patients(ext_dir) + volumes, masks, codes = [], [], [] + for p in tqdm(patients_meta, desc="Loading CROP EXT"): + try: + vol, seg, meta = load_patient_data( + p['dicom_dir'], p['nrrd_path'], verbose=False + ) + except Exception as e: + print(f" [skip] {p['code']}: {e}") + continue + if not meta.get('alignment_success', False): + continue + volumes.append(vol) + masks.append(seg) + codes.append(p['code']) + return volumes, masks, codes + + +# --------------------------------------------------------------------------- +# LOO cross-validation +# --------------------------------------------------------------------------- + +def loo_evaluate(X: np.ndarray, y: np.ndarray, task: str) -> dict: + if len(X) < 4: + print(f" [skip] {task}: only {len(X)} samples with matched labels") + return {} + + loo = LeaveOneOut() + y_true, y_pred, y_prob = [], [], [] + + for train_idx, test_idx in loo.split(X): + clf = _make_clf() + clf.fit(X[train_idx], y[train_idx]) + y_true.append(y[test_idx[0]]) + y_pred.append(clf.predict(X[test_idx])[0]) + y_prob.append(clf.predict_proba(X[test_idx])[0, 1]) + + y_true = np.array(y_true) + y_pred = np.array(y_pred) + y_prob = np.array(y_prob) + + auc = roc_auc_score(y_true, y_prob) if len(np.unique(y_true)) > 1 else float('nan') + return { + "n": len(y_true), + "accuracy": accuracy_score(y_true, y_pred), + "balanced_accuracy": balanced_accuracy_score(y_true, y_pred), + "auc": auc, + "report": classification_report(y_true, y_pred, zero_division=0), + } + + +# --------------------------------------------------------------------------- +# External validation +# --------------------------------------------------------------------------- + +def run_ext_validation( + ext_dir: str, + clfs: dict, + labels: dict, + output_dir: str, +) -> None: + print(f"\n{'=' * 60}") + print("External Validation (CROP EXT)") + print(f"{'=' * 60}") + + volumes, masks, codes = load_ext_data(ext_dir) + print(f"Loaded {len(volumes)} EXT patients") + + if len(volumes) == 0: + print("No EXT patients loaded — check ext_dir structure.") + return + + # Run predictions + rows = [] + for code, vol, msk in zip(codes, volumes, masks): + preds = predict_classifications(vol, msk, clfs) + row = {"patient_code": code} + row.update(preds) + # Attach ground-truth labels if available in Excel + if code in labels: + for task in TASK_TO_FEATURES: + row[f"{task}_gt"] = labels[code][task] + rows.append(row) + + # Save predictions CSV + os.makedirs(output_dir, exist_ok=True) + csv_path = os.path.join(output_dir, "ext_predictions.csv") + fieldnames = ["patient_code"] + list(TASK_TO_FEATURES.keys()) + \ + [f"{t}_gt" for t in TASK_TO_FEATURES] + with open(csv_path, "w", newline="") as f: + writer = csv.DictWriter(f, fieldnames=fieldnames, extrasaction="ignore") + writer.writeheader() + writer.writerows(rows) + print(f"Predictions saved → {csv_path}") + + # Compute metrics for patients that have Excel labels + labeled_rows = [r for r in rows if f"{list(TASK_TO_FEATURES)[0]}_gt" in r] + if not labeled_rows: + print("No EXT patients found in Excel — predictions saved, no metrics computed.") + return + + print(f"\nMetrics on {len(labeled_rows)} EXT patients with Excel labels:") + for task in TASK_TO_FEATURES: + y_true = np.array([r[f"{task}_gt"] for r in labeled_rows]) + y_pred = np.array([r[task] for r in labeled_rows]) + if len(np.unique(y_true)) < 2: + print(f"\n--- {task}: only one class in EXT labels, skipping AUC ---") + continue + print(f"\n--- {task} ---") + print(f" Balanced accuracy: {balanced_accuracy_score(y_true, y_pred):.3f}") + print(f" {classification_report(y_true, y_pred, zero_division=0)}") + + +# --------------------------------------------------------------------------- +# Main +# --------------------------------------------------------------------------- + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument("--data-dir", default=ExperimentConfig.DATA_DIR) + parser.add_argument("--xlsx", default=XLSX_PATH) + parser.add_argument("--output-dir", default=None) + parser.add_argument("--ext-dir", default=None, + help="Path to CROP EXT folder for external validation") + args = parser.parse_args() + + output_dir = args.output_dir or os.path.join( + os.path.dirname(__file__), '..', 'outputs', 'classification' + ) + + print("=== AEA Classification Experiment ===\n") + print(f"Data dir : {args.data_dir}") + print(f"XLSX : {args.xlsx}") + print(f"Output : {output_dir}") + if args.ext_dir: + print(f"EXT dir : {args.ext_dir}") + print() + + labels = load_classification_labels(args.xlsx) + print(f"Loaded {len(labels)} patient labels from XLSX") + + print("\nLoading imaging data...") + volumes, masks, codes = load_data(args.data_dir) + print(f"Loaded {len(volumes)} valid CROP1 patients") + + matched = [c for c in codes if c in labels] + print(f"Matched to XLSX labels: {len(matched)}/{len(codes)}") + unmatched = [c for c in codes if c not in labels] + if unmatched: + print(f" Unmatched codes: {unmatched[:10]}{'...' if len(unmatched) > 10 else ''}") + print(" → Adjust _extract_patient_code() if codes look wrong") + + if len(matched) == 0: + print("\nERROR: No patients matched. Check patient code extraction.") + return + + X_dict, y_dict = build_feature_matrix(codes, volumes, masks, labels) + n = len(y_dict["aeal_roof_contact"]) + print(f"\nFeature matrix built: {n} patients with matched labels") + + # LOO cross-validation + print(f"\n{'=' * 60}") + print("Leave-One-Out Cross-Validation (CROP1)") + print(f"{'=' * 60}") + + for task, feat_key in TASK_TO_FEATURES.items(): + X = X_dict[feat_key] + y = y_dict[task] + print(f"\n--- {task} ---") + print(f" Features : {feat_key} shape={X.shape}") + print(f" Positive : {int(y.sum())}/{len(y)} ({100*y.mean():.1f}%)") + res = loo_evaluate(X, y, task) + if res: + print(f" Accuracy : {res['accuracy']:.3f}") + print(f" Balanced accuracy : {res['balanced_accuracy']:.3f}") + print(f" AUC : {res['auc']:.3f}") + print(f" {res['report']}") + + # Train final classifiers on all CROP1 data + print(f"\n{'=' * 60}") + print("Training final classifiers on all CROP1 data...") + clfs = train_classifiers(X_dict, y_dict) + save_classifiers(clfs, output_dir) + print(f"Saved {len(clfs)} classifiers to {output_dir}/") + for task in clfs: + print(f" → {task}.joblib") + + # External validation + if args.ext_dir: + run_ext_validation(args.ext_dir, clfs, labels, output_dir) + + print("\nDone.") + + +if __name__ == "__main__": + main() diff --git a/try2/shared/classification.py b/try2/shared/classification.py new file mode 100644 index 0000000..1b80fd1 --- /dev/null +++ b/try2/shared/classification.py @@ -0,0 +1,264 @@ +""" +Post-hoc classification of AEA segmentation masks. + +Given a 3D volume (H×W×D) and a 3-class segmentation mask (BG=0, AEAL=1, AEAR=2), +predicts 4 binary labels per patient: + + Roof contact (geometric — mask z-position only): + aeal_roof_contact : AEAL in contact with skull roof (1) vs distant (0) + aear_roof_contact : AEAR in contact with skull roof (1) vs distant (0) + + Ethmoid sinus status (radiological — CT intensity in mask region): + left_ethmoid_filled : left ethmoid filled/sinusitis (1) vs clear (0) + right_ethmoid_filled : right ethmoid filled/sinusitis (1) vs clear (0) +""" + +import os +from typing import Dict, List, Optional, Tuple + +import numpy as np +import joblib +import openpyxl +from sklearn.linear_model import LogisticRegression +from sklearn.preprocessing import StandardScaler +from sklearn.pipeline import Pipeline + + +def load_classification_labels(xlsx_path: str) -> Dict[str, Dict[str, int]]: + """Return dict keyed by patient_code with 4 binary label values.""" + wb = openpyxl.load_workbook(xlsx_path) + ws = wb.active + labels: Dict[str, Dict[str, int]] = {} + for row in list(ws.iter_rows(values_only=True))[1:]: + code = row[0] + if code is None: + continue + code = str(code).strip() + labels[code] = { + "aeal_roof_contact": 1 if row[4] == "X" else 0, + "aear_roof_contact": 1 if row[6] == "X" else 0, + "left_ethmoid_filled": 1 if row[9] == "X" else 0, + "right_ethmoid_filled": 1 if row[11] == "X" else 0, + } + return labels + +def extract_roof_features(volume: np.ndarray, mask: np.ndarray, class_id: int) -> np.ndarray: + """ + Features for skull-roof contact detection. + + Slices are sorted ascending by ImagePositionPatient[2] (standard DICOM), so + high z-index = superior. "Roof contact" means the anatomy abuts the skull + base / roof, detected by sampling CT intensities just *above* (z+1..z+3) the + topmost mask voxels — bone HU >> soft tissue. + + Returns zeros if class_id is absent. + """ + voxels = np.argwhere(mask == class_id) # (N, 3): rows are (h, w, d) + D = mask.shape[2] + + if len(voxels) == 0: + return np.zeros(6, dtype=np.float32) + + z_coords = voxels[:, 2].astype(np.float32) + z_max = int(z_coords.max()) + top_z_norm = z_max / (D - 1) + centroid_z_norm = float(z_coords.mean()) / (D - 1) + z_extent_norm = (z_coords.max() - z_coords.min()) / (D - 1) + log_count = float(np.log1p(len(voxels))) + + # Sample HU values in a shell of slices just above the top of the mask. + # bone ~400–1900 HU; air ≈ -1000 HU; soft tissue -100..+100 HU. + shell_hu: List[float] = [] + top_voxels = voxels[voxels[:, 2] == z_max] # (M, 3) + for dz in (1, 2, 3): + z_above = z_max + dz + if z_above >= D: + break + hu_vals = volume[top_voxels[:, 0], top_voxels[:, 1], z_above].astype(np.float32) + shell_hu.extend(hu_vals.tolist()) + + if shell_hu: + shell_arr = np.array(shell_hu, dtype=np.float32) + bone_contact_ratio = float(np.mean(shell_arr > 400)) + mean_above_hu = float(np.mean(shell_arr)) + else: + bone_contact_ratio = 0.0 + mean_above_hu = 0.0 + + return np.array( + [top_z_norm, centroid_z_norm, z_extent_norm, log_count, + bone_contact_ratio, mean_above_hu], + dtype=np.float32, + ) + + +def extract_ethmoid_features( + volume: np.ndarray, mask: np.ndarray, class_id: int +) -> np.ndarray: + """ + 6 CT intensity features sampled from voxels belonging to class_id. + volume shape: (H, W, D), dtype typically int16 (HU values). + Returns zeros if class_id is absent. + """ + voxels = np.argwhere(mask == class_id) + + if len(voxels) == 0: + return np.zeros(6, dtype=np.float32) + + intensities = volume[voxels[:, 0], voxels[:, 1], voxels[:, 2]].astype(np.float32) + + mean_hu = float(np.mean(intensities)) + median_hu = float(np.median(intensities)) + std_hu = float(np.std(intensities)) + p10_hu = float(np.percentile(intensities, 10)) + p90_hu = float(np.percentile(intensities, 90)) + # fraction of voxels above air threshold — fluid/mucus >> -500 HU + filled_ratio = float(np.mean(intensities > -500)) + + return np.array( + [mean_hu, median_hu, std_hu, p10_hu, p90_hu, filled_ratio], + dtype=np.float32, + ) + + +def extract_features( + volume: np.ndarray, + mask: np.ndarray, +) -> Dict[str, np.ndarray]: + """Return feature vectors for both sides, for both tasks.""" + return { + "roof_L": extract_roof_features(volume, mask, class_id=1), + "roof_R": extract_roof_features(volume, mask, class_id=2), + "ethmoid_L": extract_ethmoid_features(volume, mask, class_id=1), + "ethmoid_R": extract_ethmoid_features(volume, mask, class_id=2), + } + + +def build_feature_matrix( + patient_codes: List[str], + volumes: List[np.ndarray], + masks: List[np.ndarray], + labels: Dict[str, Dict[str, int]], +) -> Tuple[Dict[str, np.ndarray], Dict[str, np.ndarray]]: + """ + Returns: + X_dict — {"roof_L": (N,4), "roof_R": (N,4), "ethmoid_L": (N,6), "ethmoid_R": (N,6)} + y_dict — {"aeal_roof_contact": (N,), "aear_roof_contact": (N,), + "left_ethmoid_filled": (N,), "right_ethmoid_filled": (N,)} + Only patients whose code appears in `labels` are included; returns their + filtered indices as well. + """ + feat_roof_L, feat_roof_R = [], [] + feat_eth_L, feat_eth_R = [], [] + y_aeal_roof, y_aear_roof = [], [] + y_left_eth, y_right_eth = [], [] + + for code, vol, msk in zip(patient_codes, volumes, masks): + if code not in labels: + continue + feats = extract_features(vol, msk) + feat_roof_L.append(feats["roof_L"]) + feat_roof_R.append(feats["roof_R"]) + feat_eth_L.append(feats["ethmoid_L"]) + feat_eth_R.append(feats["ethmoid_R"]) + lbl = labels[code] + y_aeal_roof.append(lbl["aeal_roof_contact"]) + y_aear_roof.append(lbl["aear_roof_contact"]) + y_left_eth.append(lbl["left_ethmoid_filled"]) + y_right_eth.append(lbl["right_ethmoid_filled"]) + + X_dict = { + "roof_L": np.stack(feat_roof_L) if feat_roof_L else np.empty((0, 6)), + "roof_R": np.stack(feat_roof_R) if feat_roof_R else np.empty((0, 6)), + "ethmoid_L": np.stack(feat_eth_L) if feat_eth_L else np.empty((0, 6)), + "ethmoid_R": np.stack(feat_eth_R) if feat_eth_R else np.empty((0, 6)), + } + y_dict = { + "aeal_roof_contact": np.array(y_aeal_roof, dtype=np.int32), + "aear_roof_contact": np.array(y_aear_roof, dtype=np.int32), + "left_ethmoid_filled": np.array(y_left_eth, dtype=np.int32), + "right_ethmoid_filled": np.array(y_right_eth, dtype=np.int32), + } + return X_dict, y_dict + + +def _make_clf() -> Pipeline: + return Pipeline([ + ("scaler", StandardScaler()), + ("clf", LogisticRegression(class_weight="balanced", max_iter=1000, C=1.0)), + ]) + + +TASK_TO_FEATURES = { + "aeal_roof_contact": "roof_L", + "aear_roof_contact": "roof_R", + "left_ethmoid_filled": "ethmoid_L", + "right_ethmoid_filled": "ethmoid_R", +} + + +def train_classifiers( + X_dict: Dict[str, np.ndarray], + y_dict: Dict[str, np.ndarray], +) -> Dict[str, Pipeline]: + """Train one pipeline per task. Returns dict keyed by task name.""" + clfs: Dict[str, Pipeline] = {} + for task, feat_key in TASK_TO_FEATURES.items(): + X = X_dict[feat_key] + y = y_dict[task] + if len(X) == 0: + continue + clf = _make_clf() + clf.fit(X, y) + clfs[task] = clf + return clfs + + +def save_classifiers(clfs: Dict[str, Pipeline], output_dir: str) -> None: + os.makedirs(output_dir, exist_ok=True) + for task, clf in clfs.items(): + joblib.dump(clf, os.path.join(output_dir, f"{task}.joblib")) + + +def load_classifiers(model_dir: str) -> Dict[str, Pipeline]: + clfs: Dict[str, Pipeline] = {} + for task in TASK_TO_FEATURES: + path = os.path.join(model_dir, f"{task}.joblib") + if os.path.exists(path): + clfs[task] = joblib.load(path) + return clfs + + +def predict_classifications( + volume: np.ndarray, + mask: np.ndarray, + clfs: Dict[str, Pipeline], +) -> Dict[str, int]: + """ + Given a single patient's volume and predicted mask, return classification + predictions as a dict of {task_name: 0_or_1}. + """ + feats = extract_features(volume, mask) + results: Dict[str, int] = {} + for task, feat_key in TASK_TO_FEATURES.items(): + if task not in clfs: + continue + x = feats[feat_key].reshape(1, -1) + results[task] = int(clfs[task].predict(x)[0]) + return results + + +def predict_classifications_proba( + volume: np.ndarray, + mask: np.ndarray, + clfs: Dict[str, Pipeline], +) -> Dict[str, float]: + """Same as predict_classifications but returns probability of positive class.""" + feats = extract_features(volume, mask) + results: Dict[str, float] = {} + for task, feat_key in TASK_TO_FEATURES.items(): + if task not in clfs: + continue + x = feats[feat_key].reshape(1, -1) + results[task] = float(clfs[task].predict_proba(x)[0, 1]) + return results From b4c913937b0d6bd514bc2db025eb5b147d68733e Mon Sep 17 00:00:00 2001 From: Petra Hedesiu Date: Wed, 27 May 2026 08:50:21 -0400 Subject: [PATCH 2/5] AEA classification try2 --- data_utils.py | 534 ++++++++++++++++++ .../classification/aeal_roof_contact.joblib | Bin 0 -> 1521 bytes .../classification/aear_roof_contact.joblib | Bin 0 -> 1521 bytes .../classification/ext_predictions.csv | 20 + .../classification/left_ethmoid_filled.joblib | Bin 0 -> 1521 bytes .../right_ethmoid_filled.joblib | Bin 0 -> 1521 bytes .../__pycache__/__init__.cpython-39.pyc | Bin 0 -> 462 bytes .../__pycache__/classification.cpython-39.pyc | Bin 0 -> 7693 bytes try2/shared/__pycache__/config.cpython-39.pyc | Bin 0 -> 3729 bytes .../shared/__pycache__/metrics.cpython-39.pyc | Bin 0 -> 5548 bytes 10 files changed, 554 insertions(+) create mode 100644 data_utils.py create mode 100644 try2/outputs/classification/aeal_roof_contact.joblib create mode 100644 try2/outputs/classification/aear_roof_contact.joblib create mode 100644 try2/outputs/classification/ext_predictions.csv create mode 100644 try2/outputs/classification/left_ethmoid_filled.joblib create mode 100644 try2/outputs/classification/right_ethmoid_filled.joblib create mode 100644 try2/shared/__pycache__/__init__.cpython-39.pyc create mode 100644 try2/shared/__pycache__/classification.cpython-39.pyc create mode 100644 try2/shared/__pycache__/config.cpython-39.pyc create mode 100644 try2/shared/__pycache__/metrics.cpython-39.pyc diff --git a/data_utils.py b/data_utils.py new file mode 100644 index 0000000..be7e991 --- /dev/null +++ b/data_utils.py @@ -0,0 +1,534 @@ +#!/usr/bin/env python3 +""" +Data Utilities for Medical Image Segmentation +============================================= +Handles loading and aligning DICOM images with NRRD segmentation files, +accounting for different coordinate systems used by 3D Slicer. +""" + +import numpy as np +import pydicom +import nrrd +from pathlib import Path +from typing import Tuple, Optional, Dict, List +from tqdm import tqdm +import warnings + + +def load_dicom_series(dicom_dir: str, verbose: bool = True) -> Tuple[np.ndarray, Dict]: + """ + Load a DICOM series from a directory. + + Returns: + volume: 3D numpy array (H, W, D) + metadata: Dictionary with spatial information + """ + dicom_dir = Path(dicom_dir) + dcm_files = list(dicom_dir.glob("*.dcm")) + list(dicom_dir.glob("*.DCM")) + + if not dcm_files: + raise FileNotFoundError(f"No DICOM files found in {dicom_dir}") + + if verbose: + print(f"Found {len(dcm_files)} DICOM files") + + # Load all slices + slices = [] + for f in tqdm(dcm_files, desc="Loading DICOMs", disable=not verbose): + try: + ds = pydicom.dcmread(str(f)) + slices.append(ds) + except Exception as e: + if verbose: + print(f" Warning: Failed to load {f.name}: {e}") + + if not slices: + raise ValueError("No valid DICOM files could be loaded") + + # Sort by slice position + try: + slices.sort(key=lambda s: float(s.ImagePositionPatient[2])) + sorted_by = "ImagePositionPatient" + except AttributeError: + try: + slices.sort(key=lambda s: float(s.SliceLocation)) + sorted_by = "SliceLocation" + except AttributeError: + slices.sort(key=lambda s: float(s.InstanceNumber)) + sorted_by = "InstanceNumber" + + # Stack into 3D volume, converting stored pixel values to Hounsfield + # units via the per-slice RescaleSlope / RescaleIntercept tags. HU is + # required for fixed CT windowing; min-max-normalised pipelines are + # unaffected since min-max is invariant to a positive affine rescale. + volume = np.stack([ + s.pixel_array.astype(np.float32) * float(getattr(s, "RescaleSlope", 1.0)) + + float(getattr(s, "RescaleIntercept", 0.0)) + for s in slices + ], axis=-1) + + # Extract metadata + ds = slices[0] + pixel_spacing = list(map(float, getattr(ds, "PixelSpacing", [1.0, 1.0]))) + slice_thickness = float(getattr(ds, "SliceThickness", 1.0)) + + # Calculate slice spacing from positions if available + if len(slices) > 1 and hasattr(slices[0], 'ImagePositionPatient'): + pos0 = np.array(slices[0].ImagePositionPatient) + pos1 = np.array(slices[1].ImagePositionPatient) + slice_spacing = np.linalg.norm(pos1 - pos0) + else: + slice_spacing = slice_thickness + + metadata = { + "pixel_spacing": pixel_spacing, + "slice_thickness": slice_thickness, + "slice_spacing": slice_spacing, + "rows": ds.Rows, + "columns": ds.Columns, + "num_slices": len(slices), + "modality": getattr(ds, "Modality", "Unknown"), + "patient_id": getattr(ds, "PatientID", "Unknown"), + "sorted_by": sorted_by, + "image_position": list(map(float, getattr(ds, "ImagePositionPatient", [0, 0, 0]))), + "image_orientation": list(map(float, getattr(ds, "ImageOrientationPatient", [1, 0, 0, 0, 1, 0]))), + "rescale_slope": float(getattr(ds, "RescaleSlope", 1.0)), + "rescale_intercept": float(getattr(ds, "RescaleIntercept", 0.0)), + } + + if verbose: + print(f"DICOM volume shape: {volume.shape}") + print(f"Pixel spacing: {pixel_spacing}, Slice spacing: {slice_spacing:.2f}") + + return volume, metadata + + +def load_nrrd_segmentation(nrrd_path: str, verbose: bool = True) -> Tuple[np.ndarray, Dict]: + """ + Load NRRD segmentation file with metadata. + + Returns: + segmentation: 3D numpy array with integer labels + metadata: Dictionary with spatial info and segment names + """ + if verbose: + print(f"Loading NRRD from {nrrd_path}") + + data, header = nrrd.read(nrrd_path) + + # Parse space directions (3x3 matrix) + space_directions = header.get("space directions", np.eye(3)) + if isinstance(space_directions, list): + space_directions = np.array(space_directions) + + # Parse space origin + space_origin = header.get("space origin", [0, 0, 0]) + if isinstance(space_origin, list): + space_origin = np.array(space_origin) + + # Parse segment information + segments = {} + for key, value in header.items(): + if key.startswith("Segment") and "_Name" in key: + seg_idx = key.split("_")[0].replace("Segment", "") + label_key = f"Segment{seg_idx}_LabelValue" + if label_key in header: + try: + label = int(header[label_key]) + segments[label] = value + except: + pass + + metadata = { + "space": header.get("space", "unknown"), + "space_directions": space_directions, + "space_origin": space_origin, + "segments": segments, + "header": header, + } + + if verbose: + print(f"NRRD shape: {data.shape}") + print(f"Unique labels: {np.unique(data)}") + print(f"Segments: {segments}") + + return data.astype(np.int64), metadata + + +def align_nrrd_to_dicom( + dicom_volume: np.ndarray, + nrrd_volume: np.ndarray, + nrrd_metadata: Dict, + verbose: bool = True, + dicom_metadata: Optional[Dict] = None, +) -> Tuple[np.ndarray, bool]: + """ + Align NRRD segmentation to DICOM volume using world coordinates. + + Uses DICOM spatial metadata (ImagePositionPatient, ImageOrientationPatient, + PixelSpacing, slice_spacing) and NRRD spatial metadata (space directions, + space origin) to compute a proper coordinate-based mapping. + + Returns: + aligned_nrrd: NRRD aligned to DICOM shape (H, W, D) + success: Whether alignment was successful + """ + if verbose: + print(f"\nAligning: DICOM {dicom_volume.shape} <-> NRRD {nrrd_volume.shape}") + + # Already aligned? + if dicom_volume.shape == nrrd_volume.shape: + if verbose: + print("Shapes already match!") + return nrrd_volume, True + + # --- Coordinate-based alignment (preferred) --- + if dicom_metadata is not None: + result = _align_by_coordinates( + dicom_volume, nrrd_volume, dicom_metadata, nrrd_metadata, verbose + ) + if result is not None: + return result, True + + # --- Fallback: simple permutation (shape-only) --- + result = _try_simple_alignment(dicom_volume, nrrd_volume, verbose) + if result is not None: + return result, True + + if verbose: + print("Could not find alignment") + return nrrd_volume, False + + +def _align_by_coordinates( + dicom_volume: np.ndarray, + nrrd_volume: np.ndarray, + dicom_metadata: Dict, + nrrd_metadata: Dict, + verbose: bool = True, +) -> Optional[np.ndarray]: + """Align using world-coordinate mapping between DICOM and NRRD. + + Builds the affine transforms for both volumes, computes the NRRD voxel + index corresponding to each corner of the DICOM grid, extracts the + sub-volume, and applies any necessary axis flips. + """ + # --- DICOM affine components --- + ipp = np.array(dicom_metadata.get("image_position", [0, 0, 0]), dtype=np.float64) + iop = np.array( + dicom_metadata.get("image_orientation", [1, 0, 0, 0, 1, 0]), dtype=np.float64 + ) + ps = dicom_metadata.get("pixel_spacing", [1.0, 1.0]) + row_spacing, col_spacing = float(ps[0]), float(ps[1]) + slice_spacing = float(dicom_metadata.get("slice_spacing", 1.0)) + + # DICOM ImageOrientationPatient (PS3.3 C.7.6.2.1.1): + # iop[0:3] = "direction cosines of the first row" + # = direction ALONG the row = direction of COLUMN increase + # iop[3:6] = "direction cosines of the first column" + # = direction ALONG the column = direction of ROW increase + col_dir = iop[:3] # direction of increasing column index + row_dir = iop[3:] # direction of increasing row index + slice_dir = np.cross(col_dir, row_dir) # normal to image plane + + # DICOM PixelSpacing: + # ps[0] = row spacing (distance between adjacent rows) + # ps[1] = column spacing (distance between adjacent columns) + # + # DICOM pixel (row=r, col=c, slice=s) -> world: + # world = ipp + c * col_dir * col_spacing + # + r * row_dir * row_spacing + # + s * slice_dir * slice_spacing + + # --- NRRD affine components --- + space_dirs = nrrd_metadata.get("space_directions", np.eye(3)) + if isinstance(space_dirs, list): + space_dirs = np.array(space_dirs, dtype=np.float64) + space_origin = np.array( + nrrd_metadata.get("space_origin", [0, 0, 0]), dtype=np.float64 + ) + + # NRRD voxel (i, j, k) -> world: world = space_origin + space_dirs^T @ [i, j, k] + # We invert: voxel = inv(space_dirs) @ (world - space_origin) + try: + inv_space_dirs = np.linalg.inv(space_dirs) + except np.linalg.LinAlgError: + if verbose: + print("NRRD space_directions matrix is singular, cannot invert") + return None + + H, W, D = dicom_volume.shape + + # Compute NRRD voxel index for DICOM voxel (0,0,0) and the three unit steps + world_origin = ipp # DICOM voxel (0,0,0) + nrrd_origin = inv_space_dirs @ (world_origin - space_origin) + + # Step vectors in NRRD voxel space for one DICOM voxel step + step_row = inv_space_dirs @ (row_dir * row_spacing) + step_col = inv_space_dirs @ (col_dir * col_spacing) + step_slice = inv_space_dirs @ (slice_dir * slice_spacing) + + if verbose: + print(f" NRRD origin (for DICOM 0,0,0): {nrrd_origin}") + print(f" NRRD step per DICOM row: {step_row}") + print(f" NRRD step per DICOM col: {step_col}") + print(f" NRRD step per DICOM slice: {step_slice}") + + # Build the full NRRD index array for every DICOM voxel would be huge. + # Instead, since DICOM->NRRD is an affine mapping of axis-aligned grids, + # each DICOM axis maps to exactly one NRRD axis (with possible sign flip). + # Detect which NRRD axis each DICOM axis maps to. + + steps = np.array([step_row, step_col, step_slice]) # (3, 3) + # For each DICOM axis, find the dominant NRRD axis + axis_map = {} # dicom_axis -> nrrd_axis + axis_sign = {} # dicom_axis -> +1 or -1 + for d_ax in range(3): + abs_step = np.abs(steps[d_ax]) + n_ax = int(np.argmax(abs_step)) + # Verify this is truly axis-aligned (dominant component >> others) + if abs_step[n_ax] < 1e-6: + if verbose: + print(f" DICOM axis {d_ax} has zero step in NRRD space") + return None + off_axis = np.delete(abs_step, n_ax) + if np.any(off_axis > 0.1 * abs_step[n_ax]): + if verbose: + print(f" DICOM axis {d_ax} is not axis-aligned in NRRD space: {steps[d_ax]}") + return None + axis_map[d_ax] = n_ax + axis_sign[d_ax] = 1 if steps[d_ax][n_ax] > 0 else -1 + + # Check we have a valid 1-to-1 mapping + if len(set(axis_map.values())) != 3: + if verbose: + print(f" Axis mapping is not 1-to-1: {axis_map}") + return None + + if verbose: + labels = ["row", "col", "slice"] + for d_ax in range(3): + sign = "+" if axis_sign[d_ax] > 0 else "-" + print(f" DICOM {labels[d_ax]} -> NRRD axis {axis_map[d_ax]} ({sign})") + + # Compute the NRRD index range for each DICOM axis. + # If the DICOM volume extends slightly beyond the NRRD canvas, clamp to + # the valid range and zero-pad afterwards (those edge voxels are background). + dicom_sizes = [H, W, D] + slices_per_axis = [None, None, None] # slice objects for NRRD extraction + pad_before = [0, 0, 0] # padding needed before the extracted region (per NRRD axis) + pad_after = [0, 0, 0] # padding needed after + + for d_ax in range(3): + n_ax = axis_map[d_ax] + start_nrrd = nrrd_origin[n_ax] + step = steps[d_ax][n_ax] + end_nrrd = start_nrrd + step * (dicom_sizes[d_ax] - 1) + + lo = min(start_nrrd, end_nrrd) + hi = max(start_nrrd, end_nrrd) + lo_int = int(round(lo)) + hi_int = int(round(hi)) + expected_size = hi_int - lo_int + 1 + + # Clamp to valid NRRD range, track how much padding is needed + clamped_lo = max(0, lo_int) + clamped_hi = min(nrrd_volume.shape[n_ax] - 1, hi_int) + + pb = clamped_lo - lo_int # voxels clipped at the low end + pa = hi_int - clamped_hi # voxels clipped at the high end + + # Reject if more than 30% of the axis is out of bounds + if pb + pa > 0.3 * expected_size: + if verbose: + print( + f" NRRD axis {n_ax}: range [{lo_int}, {hi_int}] has " + f"{pb + pa}/{expected_size} voxels out of bounds " + f"[0, {nrrd_volume.shape[n_ax] - 1}] (>30%)" + ) + return None + + pad_before[n_ax] = pb + pad_after[n_ax] = pa + slices_per_axis[n_ax] = slice(clamped_lo, clamped_hi + 1) + + # Extract the sub-volume from NRRD + extracted = nrrd_volume[slices_per_axis[0], slices_per_axis[1], slices_per_axis[2]] + + # Zero-pad if any edges were clipped + if any(p > 0 for p in pad_before) or any(p > 0 for p in pad_after): + pad_widths = [(pad_before[ax], pad_after[ax]) for ax in range(3)] + extracted = np.pad(extracted, pad_widths, mode="constant", constant_values=0) + if verbose: + print(f" Padded {pad_widths} to compensate for edge clipping") + + if verbose: + print(f" Extracted NRRD region: {[str(s) for s in slices_per_axis]}, shape {extracted.shape}") + + # Permute axes: we need NRRD axes in the order [axis_map[0], axis_map[1], axis_map[2]] + # so that the result is (H, W, D) matching DICOM + perm = [axis_map[0], axis_map[1], axis_map[2]] + if perm != [0, 1, 2]: + extracted = np.transpose(extracted, perm) + + # Apply flips for negative axis signs + for d_ax in range(3): + if axis_sign[d_ax] < 0: + extracted = np.flip(extracted, axis=d_ax) + + # Verify shape + if extracted.shape != dicom_volume.shape: + if verbose: + print( + f" Shape mismatch after extraction: {extracted.shape} vs {dicom_volume.shape}" + ) + return None + + # Make contiguous copy (np.flip returns a view) + extracted = np.ascontiguousarray(extracted) + + if verbose: + n_labels = np.sum(extracted > 0) + print(f" Coordinate-based alignment successful, {n_labels} label voxels") + + return extracted + + +def _try_simple_alignment( + dicom_volume: np.ndarray, + nrrd_volume: np.ndarray, + verbose: bool = True, +) -> Optional[np.ndarray]: + """Try simple axis permutations to align volumes (shape-only fallback).""" + + dicom_shape = dicom_volume.shape + nrrd_shape = nrrd_volume.shape + + permutations = [ + (0, 1, 2), + (1, 0, 2), + (0, 2, 1), + (2, 1, 0), + (1, 2, 0), + (2, 0, 1), + ] + + for perm in permutations: + transformed = np.transpose(nrrd_volume, perm) + if transformed.shape == dicom_shape: + if verbose: + print(f" Fallback alignment: transpose{perm}") + return transformed + + return None + + +def get_labeled_slice_indices(segmentation: np.ndarray) -> List[int]: + """Get indices of slices containing labels.""" + labeled = [] + for i in range(segmentation.shape[-1]): + if np.any(segmentation[:, :, i] > 0): + labeled.append(i) + return labeled + + +def load_patient_data( + dicom_dir: str, + nrrd_path: str, + verbose: bool = True +) -> Tuple[np.ndarray, np.ndarray, Dict]: + """ + Load and align a patient's DICOM and NRRD data. + + Returns: + dicom_volume: 3D DICOM volume + aligned_segmentation: Aligned NRRD segmentation + metadata: Combined metadata + """ + # Load both + dicom_vol, dicom_meta = load_dicom_series(dicom_dir, verbose) + nrrd_vol, nrrd_meta = load_nrrd_segmentation(nrrd_path, verbose) + + # Align + aligned_seg, success = align_nrrd_to_dicom( + dicom_vol, nrrd_vol, nrrd_meta, verbose, dicom_metadata=dicom_meta + ) + + if not success: + warnings.warn(f"Alignment failed for {dicom_dir}") + + # Combined metadata + metadata = { + **dicom_meta, + "segments": nrrd_meta.get("segments", {}), + "alignment_success": success, + } + + return dicom_vol, aligned_seg, metadata + + +def discover_patients(base_dir: str) -> List[Dict]: + """ + Discover all patient folders in a dataset directory. + + Returns list of dicts with 'dicom_dir' and 'nrrd_path' keys. + """ + base = Path(base_dir) + patients = [] + + for patient_dir in sorted(base.iterdir()): + if not patient_dir.is_dir() or patient_dir.name.startswith('.'): + continue + + # Find NRRD file + nrrd_files = list(patient_dir.glob("*.nrrd")) + if not nrrd_files: + continue + + # Find DICOM subdirectory — pick the one with the most .dcm files + # (e.g. patient 001 has NL001/ with DICOMs and NL001_previews/ without) + dicom_dirs = [d for d in patient_dir.iterdir() if d.is_dir()] + if not dicom_dirs: + continue + + best_dicom_dir = None + best_dcm_count = 0 + for d in dicom_dirs: + n = len(list(d.glob("*.dcm")) + list(d.glob("*.DCM"))) + if n > best_dcm_count: + best_dcm_count = n + best_dicom_dir = d + + if best_dicom_dir is None or best_dcm_count == 0: + continue + dicom_dir = best_dicom_dir + dcm_files = list(dicom_dir.glob("*.dcm")) + list(dicom_dir.glob("*.DCM")) + + patients.append({ + "patient_id": patient_dir.name, + "dicom_dir": str(dicom_dir), + "nrrd_path": str(nrrd_files[0]), + }) + + return patients + + +if __name__ == "__main__": + # Test with sample patient + import sys + + if len(sys.argv) >= 3: + dicom_dir = sys.argv[1] + nrrd_path = sys.argv[2] + + dicom, seg, meta = load_patient_data(dicom_dir, nrrd_path) + + print(f"\n=== Results ===") + print(f"DICOM shape: {dicom.shape}") + print(f"Segmentation shape: {seg.shape}") + print(f"Alignment success: {meta['alignment_success']}") + print(f"Labeled slices: {len(get_labeled_slice_indices(seg))}") + else: + print("Usage: python data_utils.py ") + + diff --git a/try2/outputs/classification/aeal_roof_contact.joblib b/try2/outputs/classification/aeal_roof_contact.joblib new file mode 100644 index 0000000000000000000000000000000000000000..3bb2772d9105d41839bfd90e2e5458afa7c6307b GIT binary patch literal 1521 zcma)5UuYaf9KJuf%k>)5L`f?qY9$y_n;c1%qEOkyJQQ3#q){T&FwXAH-LB65>C9Y` zLqb!eF_C#tGbe~a1hIKhs1_6-#46MW(V`eCD8@gOCN*MfpG1)O+udYO8y}p9-I-zL z_kHvIzHg?fekqv?;m@!;a5j)J(f|8l=2>eT; z=96G2W#6mj=blqXRIi?v$wR7GEIzmaBmC$Ip^TB-$w-)lnuT*R6A%Two1XFvY8!>H zCh~0)R8z9O7_b-)pvM^YeUK%7XO345lVZ-5)ghY$DFT84FK%Jrqi0XT++(Fc%j8JR z)7m3b<$v(BI=8cwc}%G{i5v&}pjc|?wUHfs4f3eeD(Ne?L*x!Kfm{dvG29~@t(HjXD5d0cWr7SnbC~019R5h0b)x7ZDRx7T7CPmqfEQ+(${OlpB&TnE-my- zTyE38s=V~!%H?y~z}V7wHvO@-B=79n^X|2wdfygaWVZi>m~SAyx6*yB*!_-nCf&36 zz;~y$tNO)1yZDUu=ly zGFUiR7*thbR>`zU%Szt1W&9v6poim?2rVG$RctP7Q^^}2+MFIC6>xD_C-z0NPh3m| zcr`^2DA`o1Vyfn4T6k0))^0iY^Lp1vztwkGyBdigZz$UWf*=Spi7%kt?ePKyS%lCj zJtVst6hJC_E3IMvZcX!&wJUAsYbH!)s}| zGt3i;yeso12sbg|CRA#MYlg?cya}0~#4N1I1klgA2ol1fdX5evhK22^qSXf3&4q;^ zx!2n!725_JZ{O2M=6^XpHdY$i90(J1tll3MeX9o!pb*w|y{edBg+$Es7yzD_MY}c=f SpuG-pJXQnwN9Uu1g}(v%6lw4?DBN z%Y0F%=z;OWjwoxQT3n9l0RV;Sb{mtebV$tHRvw zpd!J%>=TA82g!U&bU27%+oQ0<2!TD(?UZe=s_3py6uV5hg?%4aWyB8q@`4n}@CeKX zuX-fdX|$e{iwkGuF`3IZiNs-9C=_1Uf)Rf54P0tt>A6Ul200JsL^2=>cn_R#HEL-} zSQB}c4ysAfRR~y&2GB42*z-V^_?>B1)^v&)lV-PQ38V-J2E3Swflod)2W!_i0xc6` zF;DYB(OCQkPcyrlrO04P-3qd8?15t0h*gMyFttjKO&pH2Dkn_lJIh*o@W}dD+dsIj zz4VI;^XA3hHq^Jz9{+Usop;qYMCi@JJh)&E2Z(KTw5_e&4d=*gH&ee}c=fB%v%Z>T6TN%h z%c!4Uxpe1YdPBYIUj6pU+86546WQypPyRJh>)*kP%=W(!Gj+ttbKQSVS69?OM$eDu zo?1~K9bQ^~S~;y6@1OecW93J6;`@jCzUjNAzPbGSnY}%Hq;@$@1KKtV#ki$QasQF+pxjRui;25GLqYZ6K1U0rC~{--VHgF{a}RiLH)*-A(E)Dek_R zPtEVPbuu*B`3l;+xbLPqadl?pn!UXG+$$T&r5~=ZRv$U~)@*G}?YQ;j&)@Z|uU<;r z{GexXQ>_i|@XXX(cds0Es`V75B2WL+`dO(sw^XlV9yXs13`s->hvu@P3$3smugrt? Q+QfGKD#$;z7#&pp222`g8~^|S literal 0 HcmV?d00001 diff --git a/try2/outputs/classification/ext_predictions.csv b/try2/outputs/classification/ext_predictions.csv new file mode 100644 index 0000000..dc5c780 --- /dev/null +++ b/try2/outputs/classification/ext_predictions.csv @@ -0,0 +1,20 @@ +patient_code,aeal_roof_contact,aear_roof_contact,left_ethmoid_filled,right_ethmoid_filled,aeal_roof_contact_gt,aear_roof_contact_gt,left_ethmoid_filled_gt,right_ethmoid_filled_gt +EXT01,1,0,1,1,,,, +EXT03,0,0,0,1,,,, +EXT04,0,1,0,1,,,, +EXT05,0,0,1,1,,,, +EXT06,1,1,1,1,,,, +EXT07,0,1,0,1,,,, +EXT08,0,1,0,0,,,, +EXT09,0,0,1,0,,,, +EXT10,0,0,0,1,,,, +EXT11,0,1,1,1,,,, +EXT12,0,0,0,1,,,, +EXT13,1,0,0,1,,,, +EXT14,0,1,0,1,,,, +EXT15,1,1,1,1,,,, +EXT17,0,0,0,1,,,, +EXT18,1,1,0,1,,,, +EXT19,1,1,1,0,,,, +EXT24,1,1,0,1,,,, +EXT25,0,0,0,1,,,, diff --git a/try2/outputs/classification/left_ethmoid_filled.joblib b/try2/outputs/classification/left_ethmoid_filled.joblib new file mode 100644 index 0000000000000000000000000000000000000000..a6f7f31e8cf2bb30c60f7a981e4b72e90d3e8f9a GIT binary patch literal 1521 zcma)5UuYaf9M0vE%X!VeDz#!%EQ%yGIg_-O5~{bFhk~c4G+>AthS}ZO+ihn5I5XFH zY9tYQO<-Oma~KeWrncaNK8X0xB8k3;e-J9PP^%WJ1gj4bD?W%2f4iIPY2$`87i5-{H&oI74Qk4!AElAU7$dnRV)SbRJwb)N%owh!m>p(nQOVZm2>eY^ z4M@1t>>N}J^C#31)mu!<)FD+Wl^)xG5q*kHLg}N~9AOKmCcV>A-Hz?*@*}G(0C`CvxW~?#1lxb3c^*fyP|i7k8F%PXRi0Z zy4}Ezde6?Vgs zr0jcQTiEy6y>?&t5c__+-yX0B?csTQUx?V=K-<{p?S1l1cJPWu`(L_teeQ~udhqo5 zAI`k4?R~}j;?|e1XzLfhxqM>zs51$l$krp?)e!1KKn-Iz||Gd=J_( zls}LkQWd?XWZIyXl6P$xKa30L;rb?_c|`rP!-ZogdGq6&(<7t|E{^IXzGw@Ghp7Os zHq&7xTg(ckYC)!jPu1bV?FX*@a>wX^-G8|7Yb=87up9wF5QG^d5YX+{A#JP^k^B={^VZHp$#HW>HNhfPTh9kPr^lvvdeCEF52z>~6?zHYx(^LM}ul3MVZ*s#70!8prApCQk7Z@MW}!BWp-*=+9b~Tv@hw2oMFoL<8^{G!7b%|D)2Ue%2b>D=3T=%LLA&D>aQJAKMi5cNq z^~5(h@wsObA*frrM~A*2QKX4CVl@b7R@ ziC|gbL?G%SolmQcE)v-B8SDr`U{7s!D~?~ynjR-vr@|z~oa3sF*kP_uYt<4SfxjiI zJ_&Z3=$(3I`jkGbd-Exk+^-9T!o77E;U}Mvx;UDdjD$((DLAKE0-}J|{jgVLw3rQR zBA=R|npRzffW_Sb^bwAIA7qK&nGuzu$*>UG-mF>!DFT84FJWWglTS~=+~ZcDWpXs; zX>V0crGN0Wd)Bg4ZJ04{961j5L9w2QRfvEvqh;42u0+MEBtqvqE5(lBk&W@f@`dN! ztE0xJ=We<3MEi*G<;mQ~+de;H9DMfHF=y`3hx5OkJyn>qjLTn~tRMXA_grn)iYBW0 zIOa&(J0e@yJMCV(FL;Q(+wQjq>_L0aw0%c_*xo>^*Wdl%t%WnUeP_t!U5ly3i^hwu z(-R-{KW%LP;j8h9uP+))ug$NY?)uSq?%12Psh!`B)cRNPBD4K3#I^?F^i%t0ANs~J z1_pnfJ^!9%Jox91$I{P?7?=Bhd4KCU%Xm4l>B+==$@t_t-hIy-X0G;XoCdUFYIK}% zK|L4RGL#+84(Up9N~@H~N?L7OHGU8m(8KX8!m@~YWhx{!wc2#w%JeWPgNwsDu`gPE z;$kMjt4(Z=R?SugGd-hH(qsBS?xx|NuH8I3VD}%$Ek`1#4J9fe2$FD<_!8RP6)#Yb zMF^eFCsliBHZ)FX*_IIEB4*e%iKSILha`Fq5)wKQhKIF{7Y1-lcn zn+Xd+a&NRvQ*0Y>Jl)$!=6^XpHdYGS90(J1tTqrylz@DN{MTS4V~nG5g=D?^Uw2zf zkCfKl%%`W;+PWn)+4ubIi9fE-bsNrpmdmD>`i|XxC0Dxi&X(;< z=NDeOv}L$^Vlh`6T;TxR1_`UYK8q>)8A4y8=TJwNXNIMahF! WP}i@>puG-pJYEI)$7Z5~*}nnEEN-X( literal 0 HcmV?d00001 diff --git a/try2/shared/__pycache__/__init__.cpython-39.pyc b/try2/shared/__pycache__/__init__.cpython-39.pyc new file mode 100644 index 0000000000000000000000000000000000000000..c3c480dead0f1643027d0ce7b5ff462bdb167271 GIT binary patch literal 462 zcmaKozfQw25XSB1Z>tuSvN3hdQcGoFKnMv~VL&WitjNALMsXb338E})JO~@F~k^HbFjKj^_`r7@?0){Ed+>F8r~l9s&tmBZVo(As^Kv712n= zG*$^sR7z8o(M;tu$7qX0^oHqJ#CNDjK7DUgU^l;hYQabaO?RoAjocQgTX3z~z5^yC z2WGi8;0^}O059v>9S(Kt1!<*j+mSq6pm+yU&>6|CTj;ivN0S;W|Rm4!`7{?5HU!vRTG z#hc7vZr|tWd(S=RI}iBN(=`pB|M}I)P45*=`)_JYeyV7^g)40x!NcVI+ z&B*AQp4qiLt8069*YO-(jaNKZ&Z~M=`K)<0JgsP|Tlea6jvY;RXS^9%ccNq6hS$)w zU^Xzz-`u)!TB9GYJTkb;tN&$!m`+KE>-8Ir*HCX+DFIQ{xe$j8&56 zu}3Cv@L7z%5Lm(NX$>XDCZ5miBXH#2N(z&8NpbFx&5!dFpg12i_{m2)KgC}_>olL| zr}-KFB0u}c^3F`Id#2uylwQK#&hc+x?Td1T>{CbryVBWbuG4C642Q;tAU&B1cjXPTN^}c6lTJ57G%`R*tZPt$b zG!56owx5Mb%#wBX_RY6h8g#lroT{Y@rh*-5?WToBh|{|CG@N+1`a-l7{LKaWrUN$0El3 zz>iu2Yg?nW?7G|q3*#}l`(d`h(ye|JF`CL2E;ZRs%6OP&evG$WY>sG&=h0gIGiia? z&1|Eagq)>e+)o)^li{%-d>$rI(!pg!@?h`W1?en^ajpTqL>^QG9o$rTz9?AFT7jb0 zS`VWr;0!0EZdP_P)x9j=1qR_o>U`}e@WuFaBJ6BDL+Uz{-G2tDLjs)JO{=h0aIb}R z3rsBB4|{YmepJ|Z`@JY=>c!kGoHNbB_G-`(0Ynf-(uoWk;iuU7dXyq2SPisRvx4G5V#CNwBGt?jG`_2LNTHGNB+cL!8ax*qK z^|)nEZ`@NYd;VeEQ!9IZs(Sd;BP%y^`weY-F*l)-^}QRb@-&b*2=`Vnz$LB?M&b|2 zjEW3#*uYMl)pB zQ2{TGKS#cb6VDTsW~FeGUJ&;NyHPP6B|dN6PsG+*l57Q8upi0 zk~RP-&eCNlF3nyRSC%s|xV#+f0!l(PwzRaIZuo$l<)g5-)Eg9M86u^82{%h;L#Xa_=Z?EG$yi!j%MDAs1wt+a?hyKgaQ`;3vg( z2e&1~RnLtr10rjZ$Z{RxXhCF6h${NU)l7r%nkXTn4}s^o$!%_bhM9!+`wop@wVhRR zd*40K9_T|ScSv;SwHLKvC99&pM*ZVHSItB@b)d~_-MRSIZmR-LmBVW=B1;H5bFv7|!B1Cs$T- zH?QWkd^Vr!NIoeJHf#L{V7^57cLM-1(1@)kLK%b$(jbg%2M@9UqMXFAhGbptL}44V zAyF(%L({>Q^Fz%27f6(_QT6RW(#Cj2<6zy2k(6NH_b(?mh?7& z!7yZD%!6HagQa~E-9#)gc7D?SooC&^k5iUy1TegQmUIWqU+cq)QjaublH8h~25ec{ zv_+=&u8<}5yFFlPhfNrm5C_@prx|z9_6_FlizB(joVmh`%C9F!f} zNp{IllL1~!;()z>kG=UOP6EpVBa`-nl1U2NQnnX}BxND2IJq292?kT-ucgr1B}x4d zj)^8=OU$75gADMvd zA9YOcwd5-@HD1MQvAiE=o||q2QPkS#D+1~Tp27FL%3e!el2_U2w}AYqS5PMsyI)3KZy+pV6<{>)=MS0fT8pa6x!5`UNh=?!;KQX_(ToG0`5qupwc}uWQQJ#u>RH8UpMRGb=1Wo6-!hsQ*n`sOH^E?;xAAXP8YsD zsaKluq6TOqRS4ncUB_s1>M4y9-$Ac8HS&J}PhLezG4!FTD3NM7ASR!^OBdDiujI3? zyj~^$PYe9~NZ`kwb1CyJyi-|vcScJ5bz?@SZv&D&2f3cnt*-=GY{1bI1dz6!&};Oe zx{IH*?blGrG+022vyF$a-m(tSMja07xL?&}KerxPa1HIfYo%*gN(|e90ZN3jP-_Xr zEg?LhMQO2QJCo6xf=pdvle{*#e2$ERZk5l~qB=VLttZvhrXf{Uoza*$Bl%tvn9&%q za3sn(##5Ppa0n#ueJnmeU(+t8%Ro+gvT0FKQk)jHj4Z`DqAXcR-sOh1WUErP zreq7z0*!57%`H)6<~|~h13h+8vw8JE&u#HL*$WqLqSt}*LvyU$InrN2Ka>5H$$ma* ze}@Sf+DX?e6mnlcCPY90efW9Ge*{TbdIizqCNx02OT~Lsyidg(6{Hd31QjPyKv8is znsk-K3e^ax#VifEg}E2@r24_v!4_jKQJP{czr9u1LARF;B<>bA5cbMtI5gj*<=>~` z78RtG-V|xKatm%vHdeiw9Kj8?ZcWrzy{SQKd|_onU6EHCv__Xe17n^GhBp-^Kmm>y zLMEFiBQ}iIqQX&9Ryax`p2qPLYyB|dtrCLZYax={{U?x1{}P2ZQ`6Qq3$dB8wgTI9fTRN}5`|GAz)Rk9!KOj*oT)2%bijxcd z9atKAsUW9GP#izO1qsU1hpXcBTX6fL^6)Pz$ruYGke@w+Q zhK%Ug_(?;yiVAWl@W}S{1In4a^&?|wB5#7AQRYqb9C$f&;KLzn!pVobyo~_z0Rqpw z0{nuPbYK8!D>=M8&w%5;xR0Fo_RVk4*JXm*6gSnPL$C8-B6 zGmX1%-}!#)?uV^+Z@zu^-s;UeknJq!rphfLCKTp+m^Ixe2S#QgT2Bj$h!w6BX{2Ao zG|klqg1npwZc9eb!c`7qFk0nJk%{;Z%t}cwHQR!tKrV3|tanh>bU~MlD{SLY)RKPI z>t{eW@hq#!_)!W(vidDl;8-GYC4B~BnznC9roRHFgY6RdWhP+XA#?XH8aMydIKUSM z<`mb!SnD=YLxmgcB5-LXTfeU?Z6<3`xQ1Mfk)$MKTaZMN3LkGP95q;2yx;8=mQ1r* zn@JcqExCKKfmwo6@jSjXQh$f!Z_JioV`fV7s!9JzY@bqP0$MR< z0i;dONYXAt>W~oGHz8|uRGp3rCQu);X$NE>Z7>2gVg_W(%}D``+hP+0n-%dX9%74% zgo+*t3huy~RCK9fQ342O6*NrOrOqm6OFes{Y>gwiB8b3$VtQIfp&>iuB0@s6q>Db< zAovO%DpHhGrH+|2xBMK@F~T0@!N5Z?aegb)%Zq;l7jGk-W9J4unT73BNX5Y0*f&W; zz>P{~L9gJA?Atj!0uEiIG>^|Id5+-W_m_xcMix&5-UHyB0l>todtUME@lK+tvR;cHp6q zGA^mK$qr64av6yG$hVSr{)Fe^J~gkRfWNN(E7OvM#efE$(twvnQ$`A3vQsVj16C-e zI!`h6?=aGQxu~>SJZZODQX!FBLEIX51Vv4PG(!*&2(mkc9RrO9D#j!=6{k>mC(qae`k@V8eX_`|L--n>uEs1fte4SD!aj8b!uTCYw8R z1w8wzeh!Hgz5Zo=#(_w>NUfavL;d26Yd}2lee6_ee(lTp>=|diMswlnJg%P@)&2Y3isP1$iF2}3n0iE=L* zlI%B*s$z)7|^2^iXJ5vcs`YTq3iwu=zxV;?^mSmBJ#PCTQq)8Yd$Z~@=WhM(I4TUMT zv=iD3+mU_3qJl^4C~&OLG2nQgM(FIr|3ky~jBs;Pb-5arER%oBX zI?y8gMw_rApiK&G6#9*@QF{z5V!&}WZchLwfRk)e=tqIa*fC*04m`o81Wy7VW5wr@&n^qT1$>)b z5q<9f-(^}v~u7b5BRwD)voP= zSG4?WyS3?T@#u!>*$t;<^D*rw6|I^tY9-w$sYR`Ob_JSJ_M0B5y6DnzP2Xv^cs!@o zvekS=&zDOql&zL4w|F$C>Z*~;S9!#>y-vgTc)Yx>SJri7MXhNE7#BQVS}z*ef?BI- zwY1FW$|_S-O)Usn<#Ev@mQZtk1-{9n`QoZk%QLMzW^ZmfbqAK;vi+MIcX>i9Wy?9O zD$H5bX_{L$G{I*2`bS!Rbxp4wz&J2-zGP(A)KW<+)R^McJ)4iJx(+*rJsA4qigti; z&11Mity0de)eaP)SW)%tnpnfu&FqY^> z-l(X#T)wo*6SeG`mRm2tE7j)@D7KFcsfz0beQWIx-@W?xpSz2E;;;qgGm81r)~nxQ z`F@FyJcP~EDuuj$a4nccE#=BZqo!#&E;n~QF4uQEJe<=W=Cc}?(v-`-xy8e-?RVT( z9j6`?OjUReS-8E=K={%g*_ZaHMfS;a+9PhNN56se+9Tgczmr?EN1w~S+>>@796*L< zq@Fw@Jy&|lPS}#+Qog1{>M0Nup;9+F`*^lFYZVt#eer(JX#Hw=P@nolgpO^0B64bRP>ioz%nP`1pb%_8lN zZ$Rbul6xN-|5MDYd$#Lkw%e{_ZrhgaIh{w)+&dZHeRex@7>-_M zZZ6|(!z!&z5Seqk&-hr=d}13g$8InU&f($eBsxl_$rKqS-PGGzVD7C>Y2wa9k2-qm zpy1Z0;TC7W0W1RVTiTJ}jZ+KXHNJZ|5laDv(*O<7sY(|FJp=#S<)Y&jf1NOyyg0{vmuJE|{K2lc@-upiEWK}(>2C-L;>9m;t9 zfdczLFLK@y{UE!E94T4?`+**_-tfFf^52kskaPM&_w}5_(#J&zyAkkiH;O`pHT>X7 zZW6^PiZK+&P#i~b0>w!br%+r#aR$X%6z5Q!N0CBt69iXw+DJ#>VCc5++F+yPK(Jl zOv7aIb@#?Qd-dCpV(b3gLZc01dy8|2onda}sp)z=T+FB#QSqFAr@t+NK})@M$!fLGeEwS-AN{G`WQ0GKwoG&>P%qAb4ck^h_Ui z81mb09Zq?}tUs|^mgio_=7i%J=2O!V2Ok7RyxuT9uW9?+Z7cY!F5qa(Xs+OSx_PYQ zCv8Yzz_|E#fYIqyLBwc+Bd!rI_EVn0MFeH~LMnTk*{K0}EZ zKRxP}aW3%_0KW(RcEi~S*7*qa6%agPZg|*)Ua;xl$+MWDujOmTgM2}AKLHySH)?BY oRm+L9!v(uPMuB@6+Lzd-$Oh4cz>5}fW)IIXiu<6^!7k+g0(dcqD*ylh literal 0 HcmV?d00001 diff --git a/try2/shared/__pycache__/metrics.cpython-39.pyc b/try2/shared/__pycache__/metrics.cpython-39.pyc new file mode 100644 index 0000000000000000000000000000000000000000..6eb37ff1fb524100b1bf77c4e0f5a91745cfed04 GIT binary patch literal 5548 zcmds5O^h5z6|So8>F(+I+1c5(*K0eD6OzD$u$wp`3SevlCPa`dA$ACi5>47uy}PqJ z(><>0F}77B4h|wnSP>EzA|G~-=E4neh=hc4=8DuO5QoT@2*CyAN8)?cJLBE8vBd!n z^r&mzt9sR4_1;(Sr=waeWVn8OxUup6A;$hn%p7g=b#0}dM+z_%X zL-fo}B`UHCz9}DsOf`0J3!+xVBBNW*+;C(=&c0P@&&lKWkZH;}$jr-!?jbWT4?yOC zeAvp|?QIJ;TzNbgmBtOSmkKQ7!ue6-pJ=^PMdFlLU*pKtzyTh!PCQ-k|b5F(DNF`ku zVtNTs0q$?$*5}b=Y>RKRcOK_N3CS~OTYz@6(w1k~HowNl{52TM#qXNKS<-ll>F1t= zk?8kB+>rDQ-{$INDb6uNi*r_TueX+@D*gbkuvaU{gC9jI=_PAAuXm%YyWW}dRoxq4 z^{c&p7OTAe$|%WtgZQ;VJTI;G(dJ}b`0_ATcqNT3 zEiLQyNX2q_5NE2_)l0)Gd9Aba=8le5V*ZoB`&Nka}ahJB)IhJ`)A(VZN-Q+9KRM^wKL}yXbjq_5Ls_uI3osScps@#)t z7J+xV!%<#JMuXuMRm0(G?h;U$#mpnmHOQce1o&g=fKP4-Lb@8(bX&1vOSM z*y;8otz)f9=<56UBXp1Xv` zN7mI?Inv+NBu_Ka)1O|j)Zo4--btpdHt2#!2y|^_A`65PpTpX5E3;Vlr?Pljp~?Q#~0Zz#6=gDnCkqd z$i>o=Iq%-&c^SvDfaNZ8geUtwopsZ}@YI*^w141cS05vOl*Jv{>&6}3h3{E9*-yJs zU!MXonV=qmu=3Hg$`&Wqf_X_r$r^!AXIJpb;y_uq^XxOZqelZ(rFpdYR!z}yR-dEB zqVow*d40w|qisGg!h&klBRme;aQH?Q^*_&FVS$yMhi87Ju=FR!i@%hq|7 zTwDh{?WNfG33{77%FR1B`Y;CdQ8dg8c_3UtC%Zwd2zXr_6Nl)ee;wo9{Xis7V1f8P zi2DFBn*ebFz;TcPl3OK$-?kvY0VL%~-wZ%F?F~)Ygn$}B@lxjOfue&kg5i$DU7!e0 zVxj2nf#UBne+P=CT~PE^)E8){`=j1}fIzj32Pms@i%QbFdIrs_M0^p=)rTimP%$9Z zzc^M697c>{2%}{1o+Y`K}worr5#sOM$~)G1SUkxs~|2Sb{QAW=;u-(4B1i0IwqiaNHVyq zJ2BCY+L%bS2*pHi>WFHlrv8|DrecAv2JF6tb38hZNTH5cAdb~nTF_gZ7=H+5IXu|5d--! z=SrQX`Df7N;RKPZC-y_Jv40nrFhL(j!~7<1BB)V_gO}p-@K8jBlOXp!$aE>7O|A-v zGHyCT`uk`mh`Z(hZ;s8mEZ!_n_!9;GEbmi@+r^iGHV4@l@-GT{6tX105g-E+$lKM| z@cee&W*}3osBckcqwvMIm8HH--E-69xNGEP8mYp9(`VplQS?^USk`Q`&sJrhb=qb7 zYEo=Fuw=Kd+W}d{_W49nSJ@h_kh3H@b02TCT^K{XCRW7NMS#DD+UvT%rtzvZ#CpxpJ)6)4uH3Wb97o zu?VpMw>70fd^fISBfwSiv1{Dd9daLIKSP}Nnc`tDo z9trlF073;0p@NE6AiYMKyWu0%4mKUtM{gNvb;VS-Li|?&2Bg2^T2|SpnwqU~=+Kw9 zYNozb-$oIKr-tCF6Yk>kcht}C?Wv#LjTuwMXAbAOK=wpjxWk3AG@e1C-7o>x`Sy4g zD{qjC8#gh7N%j-eT(&J}+pfBBigRXeW45Tku)F#30W)uAH{J0;bHE(LClIl~Xg{WD z4xVIL(=@l{z=gJTq<@_b$HMrKSuiu3uKLg{jAyd4nZ?c<7e6!&06xIChY#NU6IG+IvJl1!mF`L zb*JCkjB^JfijI#uLJL?_=ZUtN9l3`hrGE_+DpJI;EO&L54ylcYc{l|r4iTM$p-M04 zc7@0y=htnoH;C5oMWQr_^d^o~I?7O(fG=d>6?N8GJ<+aPV*0?g1JZ``c7W7>=fPMk z6#~Syt2>M>)7w6}=?&@?TJ=lR*t2^3z&2BS3mWvUqHdiWtzSb-lNVt=Inj`uscEMWZ4|t|)Q0wz43xuxQWOk0MJ+va~u%Y#9~x@qsAG!`r;N zKof1IO?4ZfB<%;1;-iCu8f;fqei!#@9arwXl*$o4td&Mt4^Gl0_mzid(Y|P%n6jWD zE$mYNZ5sq9;r1~96zxZEpdl?%=xRE3zu`6ghL3M(p9t$Wn;%ynybs6O$Azrm&5x=N Xu}>}cN%b)1eN5INy9jt<@zK8k;b`## literal 0 HcmV?d00001 From 060b0b721d75df6e4b66e1f3a901781afa863b4a Mon Sep 17 00:00:00 2001 From: Petra Hedesiu Date: Wed, 27 May 2026 09:03:20 -0400 Subject: [PATCH 3/5] AEA classification try2 --- .gitignore | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..e5f29ef --- /dev/null +++ b/.gitignore @@ -0,0 +1,21 @@ +# Python +__pycache__/ +*.py[cod] +*.pyo + +# Data +CROP1/ +CROP EXT/ +*.xlsx +*.nrrd +*.dcm + +# Outputs / models +try2/outputs/ +try3/results/ +*.joblib +*.pth +*.pt + +# System +.DS_Store From 83c8070b634e79d2f6d0a0c1a022fb043c639b9e Mon Sep 17 00:00:00 2001 From: Petra Hedesiu Date: Wed, 27 May 2026 15:56:20 -0400 Subject: [PATCH 4/5] ethmoid classification on saggital plane + additional experiments --- try2/experiments/exp_classification.py | 2 +- try2/experiments/exp_ethmoid_analysis.py | 403 ++++++++++++++++++ try2/experiments/exp_ethmoid_round2.py | 361 ++++++++++++++++ .../classification/ext_predictions.csv | 20 +- .../classification/left_ethmoid_filled.joblib | Bin 1521 -> 1521 bytes .../right_ethmoid_filled.joblib | Bin 1521 -> 1521 bytes .../__pycache__/classification.cpython-39.pyc | Bin 7693 -> 8537 bytes try2/shared/classification.py | 53 ++- 8 files changed, 815 insertions(+), 24 deletions(-) create mode 100644 try2/experiments/exp_ethmoid_analysis.py create mode 100644 try2/experiments/exp_ethmoid_round2.py diff --git a/try2/experiments/exp_classification.py b/try2/experiments/exp_classification.py index 84ad008..69d2ce4 100644 --- a/try2/experiments/exp_classification.py +++ b/try2/experiments/exp_classification.py @@ -168,7 +168,7 @@ def loo_evaluate(X: np.ndarray, y: np.ndarray, task: str) -> dict: y_true, y_pred, y_prob = [], [], [] for train_idx, test_idx in loo.split(X): - clf = _make_clf() + clf = _make_clf(task) clf.fit(X[train_idx], y[train_idx]) y_true.append(y[test_idx[0]]) y_pred.append(clf.predict(X[test_idx])[0]) diff --git a/try2/experiments/exp_ethmoid_analysis.py b/try2/experiments/exp_ethmoid_analysis.py new file mode 100644 index 0000000..629e5e0 --- /dev/null +++ b/try2/experiments/exp_ethmoid_analysis.py @@ -0,0 +1,403 @@ +""" +Systematic analysis of ethmoid sinus filling prediction. +Tests HU thresholds, slab widths, distribution features, sagittal position, +classifiers, and combinations via LOO-CV. +""" + +import sys +import os + +sys.path.insert(0, os.path.join(os.path.dirname(__file__), '..', '..')) +sys.path.insert(0, os.path.join(os.path.dirname(__file__), '..')) + +import numpy as np +from tqdm import tqdm +from scipy.stats import kurtosis as scipy_kurtosis +from sklearn.model_selection import LeaveOneOut +from sklearn.metrics import roc_auc_score, balanced_accuracy_score +from sklearn.linear_model import LogisticRegression +from sklearn.ensemble import RandomForestClassifier +from sklearn.preprocessing import StandardScaler +from sklearn.pipeline import Pipeline + +from exp_classification import load_data +from shared.classification import load_classification_labels +from shared.config import ExperimentConfig + +XLSX_PATH = os.path.join(os.path.dirname(__file__), '..', '..', 'AEA_classification_sinuzite.xlsx') +DATA_DIR = ExperimentConfig.DATA_DIR + + +# --------------------------------------------------------------------------- +# Feature extractors +# --------------------------------------------------------------------------- + +def get_mask_voxels(mask, class_id): + return np.argwhere(mask == class_id) + + +def get_slab_intensities(volume, mask, class_id, w_lo, w_hi): + voxels = get_mask_voxels(mask, class_id) + if len(voxels) == 0: + return None + sel = (voxels[:, 1] >= w_lo) & (voxels[:, 1] <= w_hi) + sv = voxels[sel] + if len(sv) == 0: + return None + return volume[sv[:, 0], sv[:, 1], sv[:, 2]].astype(np.float32) + + +def get_slab_bounds(voxels, volume_shape, half_width, position='center'): + """Return (w_lo, w_hi) for a slab. + position: 'center', 'anterior', 'middle', 'posterior', or 'full' + half_width: ±half_width voxels (ignored if position='full') + """ + W = volume_shape[1] + if position == 'full': + return 0, W - 1 + + w_vals = voxels[:, 1] + w_min, w_max = int(w_vals.min()), int(w_vals.max()) + w_range = w_max - w_min + third = max(1, w_range // 3) + + if position == 'anterior': + # anterior = low W index (remember H=A-P, W=L-R ... actually W is sagittal/L-R) + # "anterior" in the A-P sense is H, not W — but task says "anterior third of W range" + # so we interpret as the first third of the mask's W extent + anchor = w_min + third // 2 + elif position == 'middle': + anchor = w_min + third + third // 2 + elif position == 'posterior': + anchor = w_min + 2 * third + third // 2 + else: # center = W centroid + anchor = int(round(float(w_vals.mean()))) + + if half_width == 0: + return anchor, anchor + return max(0, anchor - half_width), min(W - 1, anchor + half_width) + + +def base_features(intensities, threshold=-500): + """6 features matching original: mean, median, std, p10, p90, filled_ratio.""" + if intensities is None or len(intensities) == 0: + return np.zeros(6, dtype=np.float32) + return np.array([ + float(np.mean(intensities)), + float(np.median(intensities)), + float(np.std(intensities)), + float(np.percentile(intensities, 10)), + float(np.percentile(intensities, 90)), + float(np.mean(intensities > threshold)), + ], dtype=np.float32) + + +def distribution_features(intensities): + """Extra distribution shape features: kurtosis, fluid_frac, soft_frac, entropy.""" + if intensities is None or len(intensities) == 0: + return np.zeros(4, dtype=np.float32) + kurt = float(scipy_kurtosis(intensities, fisher=True)) + fluid_frac = float(np.mean((intensities >= 0) & (intensities <= 100))) + soft_frac = float(np.mean((intensities >= -100) & (intensities <= 100))) + # entropy of 20-bin histogram from -1000 to 500 HU + counts, _ = np.histogram(intensities, bins=20, range=(-1000, 500)) + counts = counts.astype(np.float64) + counts += 1e-10 # avoid log(0) + probs = counts / counts.sum() + entropy = float(-np.sum(probs * np.log(probs))) + return np.array([kurt, fluid_frac, soft_frac, entropy], dtype=np.float32) + + +# --------------------------------------------------------------------------- +# LOO-CV runner +# --------------------------------------------------------------------------- + +def make_lr(): + return Pipeline([ + ('scaler', StandardScaler()), + ('clf', LogisticRegression(class_weight='balanced', max_iter=1000, C=1.0)), + ]) + + +def make_rf(): + return Pipeline([ + ('scaler', StandardScaler()), + ('clf', RandomForestClassifier(n_estimators=100, class_weight='balanced', random_state=42)), + ]) + + +def loo_cv(X, y, clf_factory=make_lr): + loo = LeaveOneOut() + y_true, y_prob = [], [] + for train_idx, test_idx in loo.split(X): + clf = clf_factory() + clf.fit(X[train_idx], y[train_idx]) + y_true.append(y[test_idx[0]]) + y_prob.append(clf.predict_proba(X[test_idx])[0, 1]) + y_true = np.array(y_true) + y_prob = np.array(y_prob) + y_pred = (y_prob >= 0.5).astype(int) + auc = roc_auc_score(y_true, y_prob) if len(np.unique(y_true)) > 1 else float('nan') + bal_acc = balanced_accuracy_score(y_true, y_pred) + return auc, bal_acc + + +# --------------------------------------------------------------------------- +# Build feature matrices for all experiments +# --------------------------------------------------------------------------- + +def build_features_for_experiment(volumes, masks, codes, labels, class_id, + half_width, position, threshold, + use_dist_features=False): + X_list, y_list = [], [] + for vol, msk, code in zip(volumes, masks, codes): + if code not in labels: + continue + voxels = get_mask_voxels(msk, class_id) + if len(voxels) == 0: + # still need to append a zero row to keep alignment - but skip patient + continue + w_lo, w_hi = get_slab_bounds(voxels, vol.shape, half_width, position) + intensities = get_slab_intensities(vol, msk, class_id, w_lo, w_hi) + feats = base_features(intensities, threshold) + if use_dist_features: + feats = np.concatenate([feats, distribution_features(intensities)]) + X_list.append(feats) + y_list.append(labels[code]) + if not X_list: + return np.empty((0, 6)), np.empty(0, dtype=int) + return np.stack(X_list), np.array(y_list, dtype=int) + + +# --------------------------------------------------------------------------- +# Main +# --------------------------------------------------------------------------- + +def main(): + print("=== Ethmoid Sinus Filling — Systematic Analysis ===\n") + print(f"Data dir : {DATA_DIR}") + print(f"XLSX : {XLSX_PATH}\n") + + # Load data + labels_all = load_classification_labels(XLSX_PATH) + + print("Loading imaging data (this may take a while)...") + volumes, masks, codes = load_data(DATA_DIR) + print(f"Loaded {len(volumes)} patients\n") + + results = [] # (experiment_name, side, auc, bal_acc, n_pos, n_total) + + sides = [ + ('Left', 1, 'left_ethmoid_filled'), + ('Right', 2, 'right_ethmoid_filled'), + ] + + def run_exp(name, half_width, position, threshold, use_dist, clf_factory=make_lr): + for side_name, class_id, label_key in sides: + # build label dict for this task + task_labels = {code: labels_all[code][label_key] + for code in codes if code in labels_all} + X, y = build_features_for_experiment( + volumes, masks, codes, task_labels, class_id, + half_width, position, threshold, use_dist + ) + if len(X) < 4: + continue + n_pos = int(y.sum()) + n_total = len(y) + auc, bal_acc = loo_cv(X, y, clf_factory) + results.append({ + 'experiment': name, + 'side': side_name, + 'auc': auc, + 'bal_acc': bal_acc, + 'n_pos': n_pos, + 'n_total': n_total, + }) + print(f" {name} | {side_name:5s} | AUC={auc:.3f} | BalAcc={bal_acc:.3f} | pos={n_pos}/{n_total}") + + # ----------------------------------------------------------------------- + # BASELINE (reproduce original) + print("\n--- BASELINE (slab=±2, threshold=-500, LR) ---") + run_exp("Baseline(slab±2,thr=-500,LR)", half_width=2, position='center', + threshold=-500, use_dist=False, clf_factory=make_lr) + + # ----------------------------------------------------------------------- + # A. HU threshold sweep + print("\n--- A. HU Threshold Sweep (slab=±2, center) ---") + for thr in [-800, -500, -200, 0, 100]: + run_exp(f"A_thr={thr}", half_width=2, position='center', + threshold=thr, use_dist=False, clf_factory=make_lr) + + # ----------------------------------------------------------------------- + # B. Sagittal slab width sweep + print("\n--- B. Slab Width Sweep (threshold=-500, center) ---") + slab_configs = [ + ('±0(single)', 0), + ('±1', 1), + ('±2', 2), + ('±5', 5), + ('full3D', 999), # we'll handle full3D via position='full' + ] + for label, hw in slab_configs: + pos = 'full' if hw == 999 else 'center' + hw_eff = 0 if hw == 999 else hw + run_exp(f"B_slab={label}", half_width=hw_eff, position=pos, + threshold=-500, use_dist=False, clf_factory=make_lr) + + # ----------------------------------------------------------------------- + # C. Distribution shape features (slab=±2, threshold=-500) + print("\n--- C. Distribution Shape Features (slab=±2, thr=-500) ---") + run_exp("C_base_only", half_width=2, position='center', + threshold=-500, use_dist=False, clf_factory=make_lr) + run_exp("C_base+dist", half_width=2, position='center', + threshold=-500, use_dist=True, clf_factory=make_lr) + + # ----------------------------------------------------------------------- + # D. Sagittal position (slab=±2, threshold=-500) + # Also find which position gives highest filled_ratio variance + print("\n--- D. Sagittal Position (slab=±2, thr=-500) ---") + + # Compute variance of filled_ratio across patients for each position to pick best + positions_to_test = ['anterior', 'middle', 'posterior', 'center'] + best_pos = {} + for side_name, class_id, label_key in sides: + task_labels = {code: labels_all[code][label_key] + for code in codes if code in labels_all} + best_var = -1 + best_p = 'center' + for pos in positions_to_test: + ratios = [] + for vol, msk, code in zip(volumes, masks, codes): + if code not in task_labels: + continue + voxels = get_mask_voxels(msk, class_id) + if len(voxels) == 0: + continue + w_lo, w_hi = get_slab_bounds(voxels, vol.shape, 2, pos) + intensities = get_slab_intensities(vol, msk, class_id, w_lo, w_hi) + if intensities is not None and len(intensities) > 0: + ratios.append(float(np.mean(intensities > -500))) + if ratios: + v = float(np.var(ratios)) + if v > best_var: + best_var = v + best_p = pos + best_pos[side_name] = best_p + print(f" Best position for {side_name}: {best_p} (filled_ratio variance={best_var:.4f})") + + for pos in positions_to_test: + run_exp(f"D_pos={pos}", half_width=2, position=pos, + threshold=-500, use_dist=False, clf_factory=make_lr) + + # ----------------------------------------------------------------------- + # E. Classifier comparison (slab=±2, threshold=-500, center, no dist) + print("\n--- E. Classifier Comparison (slab=±2, thr=-500, center, no dist) ---") + run_exp("E_LR", half_width=2, position='center', + threshold=-500, use_dist=False, clf_factory=make_lr) + run_exp("E_RF", half_width=2, position='center', + threshold=-500, use_dist=False, clf_factory=make_rf) + + # ----------------------------------------------------------------------- + # F. Best combination + # From the results so far, pick best slab width, threshold, position + # Then test both classifiers with dist features + print("\n--- F. Best Combination ---") + + # Gather AUCs from experiments B (slab width) and A (threshold) for each side + # to determine best settings programmatically + def best_setting(prefix, param_key): + subset = [r for r in results if r['experiment'].startswith(prefix)] + if not subset: + return None + best = max(subset, key=lambda r: r['auc'] if not np.isnan(r['auc']) else -1) + return best['experiment'] + + # Pick best slab from B experiments (average over sides) + b_exps = {} + for r in results: + if r['experiment'].startswith('B_'): + key = r['experiment'] + b_exps.setdefault(key, []).append(r['auc']) + b_avg = {k: np.nanmean(v) for k, v in b_exps.items()} + best_slab_exp = max(b_avg, key=lambda k: b_avg[k]) if b_avg else 'B_slab=±2' + # extract half_width from best_slab_exp name + slab_map = {'±0(single)': (0, 'center'), '±1': (1, 'center'), '±2': (2, 'center'), + '±5': (5, 'center'), 'full3D': (0, 'full')} + best_slab_label = best_slab_exp.replace('B_slab=', '') + best_hw, best_pos_f = slab_map.get(best_slab_label, (2, 'center')) + print(f" Best slab from B: {best_slab_exp} (hw={best_hw}, pos={best_pos_f})") + + # Pick best threshold from A experiments + a_exps = {} + for r in results: + if r['experiment'].startswith('A_'): + key = r['experiment'] + a_exps.setdefault(key, []).append(r['auc']) + a_avg = {k: np.nanmean(v) for k, v in a_exps.items()} + best_thr_exp = max(a_avg, key=lambda k: a_avg[k]) if a_avg else 'A_thr=-500' + best_thr = int(best_thr_exp.replace('A_thr=', '')) + print(f" Best threshold from A: {best_thr_exp} (thr={best_thr})") + + # Best position from D + d_exps = {} + for r in results: + if r['experiment'].startswith('D_'): + key = r['experiment'] + d_exps.setdefault(key, []).append(r['auc']) + d_avg = {k: np.nanmean(v) for k, v in d_exps.items()} + best_pos_exp = max(d_avg, key=lambda k: d_avg[k]) if d_avg else 'D_pos=center' + best_pos_d = best_pos_exp.replace('D_pos=', '') + print(f" Best position from D: {best_pos_exp} (pos={best_pos_d})") + + # Combination experiments + for use_dist in [False, True]: + for clf_name, clf_f in [('LR', make_lr), ('RF', make_rf)]: + dist_tag = '+dist' if use_dist else '' + name = f"F_slab={best_slab_label},thr={best_thr},pos={best_pos_d}{dist_tag},{clf_name}" + run_exp(name, half_width=best_hw, position=best_pos_d, + threshold=best_thr, use_dist=use_dist, clf_factory=clf_f) + + # ----------------------------------------------------------------------- + # Summary table + print("\n" + "=" * 100) + print("SUMMARY TABLE — All Experiments Ranked by Mean AUC (averaged over Left+Right)") + print("=" * 100) + + # Aggregate per experiment (mean over sides) + exp_summary = {} + for r in results: + key = r['experiment'] + exp_summary.setdefault(key, []).append(r) + + rows_summary = [] + for exp_name, exp_results in exp_summary.items(): + aucs = [r['auc'] for r in exp_results] + bals = [r['bal_acc'] for r in exp_results] + mean_auc = float(np.nanmean(aucs)) + mean_bal = float(np.nanmean(bals)) + # individual side info + sides_info = {r['side']: r for r in exp_results} + rows_summary.append((mean_auc, mean_bal, exp_name, sides_info)) + + rows_summary.sort(key=lambda x: x[0], reverse=True) + + header = f"{'Rank':>4} {'Experiment':<55} {'MeanAUC':>7} {'MeanBalAcc':>10} {'Left_AUC':>8} {'Right_AUC':>9} {'pos/n':>7}" + print(header) + print("-" * len(header)) + for rank, (mean_auc, mean_bal, exp_name, sides_info) in enumerate(rows_summary, 1): + left_r = sides_info.get('Left', {}) + right_r = sides_info.get('Right', {}) + left_auc = f"{left_r.get('auc', float('nan')):.3f}" if left_r else ' N/A ' + right_auc = f"{right_r.get('auc', float('nan')):.3f}" if right_r else ' N/A ' + # use left side for pos/n (representative) + n_pos = left_r.get('n_pos', '?') + n_total = left_r.get('n_total', '?') + pos_n = f"{n_pos}/{n_total}" + print(f"{rank:>4} {exp_name:<55} {mean_auc:>7.3f} {mean_bal:>10.3f} {left_auc:>8} {right_auc:>9} {pos_n:>7}") + + print("\nDone.") + + +if __name__ == '__main__': + main() diff --git a/try2/experiments/exp_ethmoid_round2.py b/try2/experiments/exp_ethmoid_round2.py new file mode 100644 index 0000000..3a08ca1 --- /dev/null +++ b/try2/experiments/exp_ethmoid_round2.py @@ -0,0 +1,361 @@ +""" +Round 2 ethmoid sinus filling prediction experiments. + +Base config: middle-W slab ±1, threshold -200 HU, LogisticRegression (balanced). + +Experiments: + G. I-S stratification — superior vs inferior half of D axis + H. Texture / distribution shape — kurtosis, skew, fluid/mucus fractions, bimodality + I. HU histogram as features — 10-bin histogram + PCA(3) before LogReg + J. Sinus air volume proxy — fractions below -800, -500, -200 HU + K. Left-right symmetry — joint classifier with cross-side symmetry feature + L. Combination — best features from G+H+J together +""" + +import sys +import os + +sys.path.insert(0, os.path.join(os.path.dirname(__file__), '..', '..')) +sys.path.insert(0, os.path.join(os.path.dirname(__file__), '..')) + +import numpy as np +from scipy import stats as scipy_stats +from sklearn.model_selection import LeaveOneOut +from sklearn.metrics import roc_auc_score, balanced_accuracy_score +from sklearn.linear_model import LogisticRegression +from sklearn.preprocessing import StandardScaler +from sklearn.pipeline import Pipeline +from sklearn.decomposition import PCA + +from exp_classification import load_data +from shared.classification import load_classification_labels +from shared.config import ExperimentConfig + +XLSX_PATH = os.path.join(os.path.dirname(__file__), '..', '..', 'AEA_classification_sinuzite.xlsx') +DATA_DIR = ExperimentConfig.DATA_DIR + +# Base config constants +BASE_HALF_WIDTH = 1 +BASE_THRESHOLD = -200 + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + +def get_slab_voxels_and_intensities(volume, mask, class_id): + """Return (slab_voxels, intensities) using base config: middle-W slab ±1.""" + voxels = np.argwhere(mask == class_id) + if len(voxels) == 0: + return None, None + w_min, w_max = int(voxels[:, 1].min()), int(voxels[:, 1].max()) + w_mid = (w_min + w_max) // 2 + w_lo = max(0, w_mid - BASE_HALF_WIDTH) + w_hi = min(volume.shape[1] - 1, w_mid + BASE_HALF_WIDTH) + sel = (voxels[:, 1] >= w_lo) & (voxels[:, 1] <= w_hi) + sv = voxels[sel] + if len(sv) == 0: + return None, None + intensities = volume[sv[:, 0], sv[:, 1], sv[:, 2]].astype(np.float32) + return sv, intensities + + +def make_lr(): + return Pipeline([ + ('scaler', StandardScaler()), + ('clf', LogisticRegression(class_weight='balanced', max_iter=1000, C=1.0)), + ]) + + +def loo_cv(X, y, clf_factory=make_lr): + loo = LeaveOneOut() + y_true, y_prob = [], [] + for train_idx, test_idx in loo.split(X): + clf = clf_factory() + clf.fit(X[train_idx], y[train_idx]) + y_true.append(y[test_idx[0]]) + y_prob.append(clf.predict_proba(X[test_idx])[0, 1]) + y_true = np.array(y_true) + y_prob = np.array(y_prob) + y_pred = (y_prob >= 0.5).astype(int) + auc = roc_auc_score(y_true, y_prob) if len(np.unique(y_true)) > 1 else float('nan') + bal_acc = balanced_accuracy_score(y_true, y_pred) + return auc, bal_acc + + +# --------------------------------------------------------------------------- +# Feature extractors +# --------------------------------------------------------------------------- + +def features_G(volume, mask, class_id): + """G. I-S stratification: superior vs inferior half of D axis + gradient.""" + sv, intensities = get_slab_voxels_and_intensities(volume, mask, class_id) + if sv is None: + return np.zeros(5, dtype=np.float32) + d_coords = sv[:, 2] + d_min, d_max = int(d_coords.min()), int(d_coords.max()) + d_mid = (d_min + d_max) // 2 + inf_mask = d_coords <= d_mid + sup_mask = d_coords > d_mid + inf_ints = intensities[inf_mask] + sup_ints = intensities[sup_mask] + filled_ratio_inf = float(np.mean(inf_ints > BASE_THRESHOLD)) if len(inf_ints) > 0 else 0.0 + filled_ratio_sup = float(np.mean(sup_ints > BASE_THRESHOLD)) if len(sup_ints) > 0 else 0.0 + filled_ratio_all = float(np.mean(intensities > BASE_THRESHOLD)) + gradient = filled_ratio_inf - filled_ratio_sup # positive = more filling inferiorly + # fraction of voxels in inferior half + inf_frac = float(len(inf_ints)) / max(1, len(intensities)) + return np.array([filled_ratio_inf, filled_ratio_sup, filled_ratio_all, gradient, inf_frac], + dtype=np.float32) + + +def features_H(volume, mask, class_id): + """H. Texture / distribution shape features.""" + _, intensities = get_slab_voxels_and_intensities(volume, mask, class_id) + if intensities is None or len(intensities) < 3: + return np.zeros(8, dtype=np.float32) + mean_hu = float(np.mean(intensities)) + std_hu = float(np.std(intensities)) + filled_ratio = float(np.mean(intensities > BASE_THRESHOLD)) + skew = float(scipy_stats.skew(intensities)) + kurt = float(scipy_stats.kurtosis(intensities, fisher=True)) + fluid_frac = float(np.mean((intensities >= 0) & (intensities <= 100))) + mucus_frac = float(np.mean((intensities >= -100) & (intensities < 0))) + # bimodality coefficient + n = len(intensities) + if kurt > -1 and n > 3: + bimodality = (skew ** 2 + 1) / (kurt + 3 * (n - 1) ** 2 / ((n - 2) * (n - 3) + 1e-9)) + else: + bimodality = 0.0 + return np.array([mean_hu, std_hu, filled_ratio, skew, kurt, fluid_frac, mucus_frac, bimodality], + dtype=np.float32) + + +def features_I_raw(volume, mask, class_id): + """I. 10-bin HU histogram features (raw, before PCA).""" + _, intensities = get_slab_voxels_and_intensities(volume, mask, class_id) + if intensities is None or len(intensities) == 0: + return np.zeros(10, dtype=np.float32) + counts, _ = np.histogram(intensities, bins=10, range=(-1000, 500)) + total = max(1, counts.sum()) + return (counts / total).astype(np.float32) + + +def features_J(volume, mask, class_id): + """J. Sinus air volume proxy: fractions below -800, -500, -200 HU.""" + _, intensities = get_slab_voxels_and_intensities(volume, mask, class_id) + if intensities is None or len(intensities) == 0: + return np.zeros(3, dtype=np.float32) + frac_air = float(np.mean(intensities < -800)) + frac_below500 = float(np.mean(intensities < -500)) + frac_below200 = float(np.mean(intensities < -200)) + return np.array([frac_air, frac_below500, frac_below200], dtype=np.float32) + + +def features_GHJ(volume, mask, class_id): + """L. Combined G + H + J features.""" + g = features_G(volume, mask, class_id) + h = features_H(volume, mask, class_id) + j = features_J(volume, mask, class_id) + return np.concatenate([g, h, j]) + + +# --------------------------------------------------------------------------- +# Build feature matrices +# --------------------------------------------------------------------------- + +def build_X_y(volumes, masks, codes, labels_all, class_id, label_key, feat_fn): + X_list, y_list = [], [] + for vol, msk, code in zip(volumes, masks, codes): + if code not in labels_all: + continue + feats = feat_fn(vol, msk, class_id) + X_list.append(feats) + y_list.append(labels_all[code][label_key]) + if not X_list: + return np.empty((0, 1)), np.empty(0, dtype=int) + return np.stack(X_list), np.array(y_list, dtype=int) + + +def build_X_y_joint(volumes, masks, codes, labels_all): + """K. Build joint feature matrix with both sides + symmetry features.""" + X_list, y_L_list, y_R_list = [], [], [] + for vol, msk, code in zip(volumes, masks, codes): + if code not in labels_all: + continue + # Base features for both sides + _, ints_L = get_slab_voxels_and_intensities(vol, msk, class_id=1) + _, ints_R = get_slab_voxels_and_intensities(vol, msk, class_id=2) + + fr_L = float(np.mean(ints_L > BASE_THRESHOLD)) if ints_L is not None and len(ints_L) > 0 else 0.0 + fr_R = float(np.mean(ints_R > BASE_THRESHOLD)) if ints_R is not None and len(ints_R) > 0 else 0.0 + + mean_L = float(np.mean(ints_L)) if ints_L is not None and len(ints_L) > 0 else 0.0 + mean_R = float(np.mean(ints_R)) if ints_R is not None and len(ints_R) > 0 else 0.0 + + # Symmetry features + ratio = fr_L / (fr_R + 1e-6) + diff = fr_L - fr_R + mean_diff = mean_L - mean_R + + feats = np.array([fr_L, fr_R, mean_L, mean_R, ratio, diff, mean_diff], dtype=np.float32) + X_list.append(feats) + y_L_list.append(labels_all[code]['left_ethmoid_filled']) + y_R_list.append(labels_all[code]['right_ethmoid_filled']) + + if not X_list: + return np.empty((0, 7)), np.empty(0, dtype=int), np.empty(0, dtype=int) + return np.stack(X_list), np.array(y_L_list, dtype=int), np.array(y_R_list, dtype=int) + + +# --------------------------------------------------------------------------- +# PCA-LR factory for experiment I +# --------------------------------------------------------------------------- + +def make_pca_lr(n_components=3): + def factory(): + return Pipeline([ + ('scaler', StandardScaler()), + ('pca', PCA(n_components=n_components)), + ('clf', LogisticRegression(class_weight='balanced', max_iter=1000, C=1.0)), + ]) + return factory + + +# --------------------------------------------------------------------------- +# Main +# --------------------------------------------------------------------------- + +def main(): + print("=== Ethmoid Sinus Filling — Round 2 Experiments ===\n") + print(f"Base config: middle-W slab ±{BASE_HALF_WIDTH}, threshold {BASE_THRESHOLD} HU, LogReg balanced\n") + print(f"Data dir : {DATA_DIR}") + print(f"XLSX : {XLSX_PATH}\n") + + labels_all = load_classification_labels(XLSX_PATH) + + print("Loading imaging data...") + volumes, masks, codes = load_data(DATA_DIR) + print(f"Loaded {len(volumes)} patients\n") + + results = [] + + def report(name, side, auc, bal_acc, n_pos, n_total): + results.append({'experiment': name, 'side': side, 'auc': auc, + 'bal_acc': bal_acc, 'n_pos': n_pos, 'n_total': n_total}) + tag = f" {name} | {side:5s} | AUC={auc:.3f} | BalAcc={bal_acc:.3f} | pos={n_pos}/{n_total}" + print(tag) + + def run_both_sides(name, feat_fn, clf_factory=make_lr): + for side_name, class_id, label_key in [ + ('Left', 1, 'left_ethmoid_filled'), + ('Right', 2, 'right_ethmoid_filled'), + ]: + X, y = build_X_y(volumes, masks, codes, labels_all, class_id, label_key, feat_fn) + if len(X) < 4: + print(f" {name} | {side_name} | SKIPPED (n={len(X)})") + continue + n_pos = int(y.sum()) + n_total = len(y) + auc, bal_acc = loo_cv(X, y, clf_factory) + report(name, side_name, auc, bal_acc, n_pos, n_total) + + # ----------------------------------------------------------------------- + # G. I-S stratification + print("\n--- G. I-S Stratification (superior vs inferior half of D axis) ---") + run_both_sides("G_IS_stratification", features_G) + + # ----------------------------------------------------------------------- + # H. Texture / distribution shape + print("\n--- H. Texture / Distribution Shape Features ---") + run_both_sides("H_texture_dist", features_H) + + # ----------------------------------------------------------------------- + # I. HU histogram + PCA + print("\n--- I. HU Histogram (10 bins) + PCA(3) + LogReg ---") + # Check n_components doesn't exceed n_features or n_samples-1 + # We have 10 features; use PCA(3) + run_both_sides("I_histogram_PCA3", features_I_raw, clf_factory=make_pca_lr(n_components=3)) + + # ----------------------------------------------------------------------- + # J. Sinus air volume proxy + print("\n--- J. Sinus Air Volume Proxy (fractions <-800, <-500, <-200 HU) ---") + run_both_sides("J_air_proxy", features_J) + + # ----------------------------------------------------------------------- + # K. Left-right symmetry (joint classifier) + print("\n--- K. Left-Right Symmetry (joint features, separate LOO-CV per side) ---") + X_joint, y_L, y_R = build_X_y_joint(volumes, masks, codes, labels_all) + if len(X_joint) >= 4: + auc_L, bal_L = loo_cv(X_joint, y_L) + report("K_LR_symmetry_joint", 'Left', auc_L, bal_L, int(y_L.sum()), len(y_L)) + auc_R, bal_R = loo_cv(X_joint, y_R) + report("K_LR_symmetry_joint", 'Right', auc_R, bal_R, int(y_R.sum()), len(y_R)) + else: + print(f" K_LR_symmetry_joint | SKIPPED (n={len(X_joint)})") + + # ----------------------------------------------------------------------- + # L. Combination G + H + J + print("\n--- L. Combination: G + H + J features ---") + run_both_sides("L_GHJ_combined", features_GHJ) + + # ----------------------------------------------------------------------- + # Summary + print("\n" + "=" * 110) + print("ROUND 2 SUMMARY — Ranked by Mean AUC (Left+Right average)") + print("=" * 110) + + exp_summary = {} + for r in results: + exp_summary.setdefault(r['experiment'], []).append(r) + + rows = [] + for exp_name, exp_results in exp_summary.items(): + aucs = [r['auc'] for r in exp_results] + bals = [r['bal_acc'] for r in exp_results] + mean_auc = float(np.nanmean(aucs)) + mean_bal = float(np.nanmean(bals)) + sides_info = {r['side']: r for r in exp_results} + rows.append((mean_auc, mean_bal, exp_name, sides_info)) + + rows.sort(key=lambda x: x[0], reverse=True) + + ROUND1_BASELINE = {'Left': 0.508, 'Right': 0.688, 'Mean': 0.598} + + header = (f"{'Rank':>4} {'Experiment':<30} {'MeanAUC':>7} {'BalAcc':>7}" + f" {'Left_AUC':>8} {'Right_AUC':>9} {'Left_vs_R1':>10} {'n':>6}") + print(header) + print("-" * len(header)) + + for rank, (mean_auc, mean_bal, exp_name, sides_info) in enumerate(rows, 1): + left_r = sides_info.get('Left', {}) + right_r = sides_info.get('Right', {}) + left_auc = left_r.get('auc', float('nan')) + right_auc = right_r.get('auc', float('nan')) + left_diff = left_auc - ROUND1_BASELINE['Left'] + diff_tag = f"{left_diff:+.3f}" + n_pos = left_r.get('n_pos', '?') + n_total = left_r.get('n_total', '?') + print(f"{rank:>4} {exp_name:<30} {mean_auc:>7.3f} {mean_bal:>7.3f}" + f" {left_auc:>8.3f} {right_auc:>9.3f} {diff_tag:>10} {n_pos}/{n_total}") + + print() + print("Round 1 best: Mean=0.598, Left=0.508, Right=0.688 (RF, slab±1, thr=-200, middle-W)") + print() + + # Highlight experiments that improved Left AUC + improved_left = [(mean_auc, exp_name, sides_info['Left'].get('auc', float('nan'))) + for mean_auc, _, exp_name, sides_info in rows + if 'Left' in sides_info and sides_info['Left'].get('auc', 0) > ROUND1_BASELINE['Left']] + if improved_left: + print("Experiments that IMPROVED Left AUC vs Round 1 (baseline=0.508):") + for mean_auc, exp_name, left_auc in sorted(improved_left, key=lambda x: -x[2]): + print(f" {exp_name:<30} Left AUC={left_auc:.3f} (+{left_auc - ROUND1_BASELINE['Left']:.3f})") + else: + print("No experiment exceeded the Round 1 Left AUC baseline of 0.508.") + + print("\nDone.") + + +if __name__ == '__main__': + main() diff --git a/try2/outputs/classification/ext_predictions.csv b/try2/outputs/classification/ext_predictions.csv index dc5c780..92cc202 100644 --- a/try2/outputs/classification/ext_predictions.csv +++ b/try2/outputs/classification/ext_predictions.csv @@ -1,20 +1,20 @@ patient_code,aeal_roof_contact,aear_roof_contact,left_ethmoid_filled,right_ethmoid_filled,aeal_roof_contact_gt,aear_roof_contact_gt,left_ethmoid_filled_gt,right_ethmoid_filled_gt -EXT01,1,0,1,1,,,, +EXT01,1,0,1,0,,,, EXT03,0,0,0,1,,,, -EXT04,0,1,0,1,,,, -EXT05,0,0,1,1,,,, +EXT04,0,1,1,1,,,, +EXT05,0,0,0,1,,,, EXT06,1,1,1,1,,,, -EXT07,0,1,0,1,,,, +EXT07,0,1,1,0,,,, EXT08,0,1,0,0,,,, -EXT09,0,0,1,0,,,, -EXT10,0,0,0,1,,,, +EXT09,0,0,0,0,,,, +EXT10,0,0,1,1,,,, EXT11,0,1,1,1,,,, EXT12,0,0,0,1,,,, -EXT13,1,0,0,1,,,, +EXT13,1,0,1,1,,,, EXT14,0,1,0,1,,,, -EXT15,1,1,1,1,,,, +EXT15,1,1,1,0,,,, EXT17,0,0,0,1,,,, -EXT18,1,1,0,1,,,, +EXT18,1,1,1,1,,,, EXT19,1,1,1,0,,,, EXT24,1,1,0,1,,,, -EXT25,0,0,0,1,,,, +EXT25,0,0,1,1,,,, diff --git a/try2/outputs/classification/left_ethmoid_filled.joblib b/try2/outputs/classification/left_ethmoid_filled.joblib index a6f7f31e8cf2bb30c60f7a981e4b72e90d3e8f9a..73b2189749a6f4e9b859a1175d370180d708a068 100644 GIT binary patch delta 234 zcmVys}6K`>X#nZ&RR;6U`>m2BB7;6SU}{iAubm_W#}4@Sq^-$3vp8eY7g z<3Py=EaIrTt3Q*V0zfZE?t$5rsb)afN7>hkzh*#3sQOK_;X^r@)#E+UBF~qAaaZ-$3uy+}nqu z;y{=w4O51??>>{D0zfZ;=7yu8S!O`T2HURTb7nvaF>?nCbVEQhq?m?!)jyjJDMBt` kyFa1TiGXXu(Z5<(qLaWO=96y)K?he~ky-rdwv)F7I-(G59smFU diff --git a/try2/outputs/classification/right_ethmoid_filled.joblib b/try2/outputs/classification/right_ethmoid_filled.joblib index 3c517bfc0f1adda2fe1a0408e6adc94a94b04969..83fbd90b5799724f484080d0eb2e65e6bc1b62ae 100644 GIT binary patch delta 234 zcmVo%pP-9Y2IzD*?a-9R6hv!xavk3f@M4@tz|&)E1m_b zPG>-Gj&d756vDF)1MvX?4wC@|KQB!wZjsn0cfao=43!Jj?Z4c;(rR!v#Xs~#$FW4` kUcV+^!O5y;*FRjO4Jm|0-;-|zK?lT7@MPiD$xVc+CZG9X@J5kjX)yZ(|3W?+CbK&tqX_| z;6RsNsg;H%*glh>0zfZ_i#eiJ*JMB~Fzuqgwq!t5`y+b>jy^!)EbO}%xMV<}01}D- zt!6;P_+cwTsAIDa1MvX?6O#c2KQG*JuYvdYq(5o@g=RE9>%WzuzX=Vr%Rep#)+t^m k<3DEKu@)lMzQ3N|7B&-s(UWflK?l(}X`(~Nt&_I}I;v4`Hvj+t diff --git a/try2/shared/__pycache__/classification.cpython-39.pyc b/try2/shared/__pycache__/classification.cpython-39.pyc index 309cf80cf119c1432b956ea2dd303f99e1d5ed49..25935418b77cc9c5364f33a653652286126ad7b2 100644 GIT binary patch delta 3348 zcmZuzOKcn06`eO6lEbek>O=kYWZAOlN{THtX_2@={Y%uutpmHRTB0dC;(L_F9L}(B zhWZ%8(x@6_kpyTkExM?JM(8T5A~muo3Kacpy6K|msn|7XZO^}(olrIP-G%iwBP?$9*Yl-3QTPCl9KC~8 zbyUHZ>H8{6u;dT5(h%404a;{Bk)`%EmS!2y6rQ5VvOds^9Ba@sjFLWc`w`a9@}O^z za*Yp+DuAHu`96^dkICeYT=de?<-mTF4YDC{IR=69Y#7omuo2+L*(e)h<7|RW?rWuq z-jO`gKFLZd^3vH9kan*;guXJ#z5(uZ78 zb{^Jjo)y84O1mP=rzZr?dPb&zh8Ni9!0-%H?k%>Oo{SQbOW!`Rse@&!; zEk@)=%m<$)-%OnYiI_t;AMB<+oI8W0gjt}%A->}YqvDl?<5bI4Zg@>dVkvkneS^LT z29*dv_sIhhhxZ!1?%O2zD8113qJD&FOnXCRQJmxZu>)m~JRyyvR8?A9E9z^0w4kuq zn9_=QaYpF_?iBYAM z^i)5&M;NK8A`JeL9%v_zQY4t7-{m@EyIc%))d)->;VJjD1JX))Nl$q{<);qF9%-c| zmp&UlTIVk{FEKvj3VA2mk`uwz?@w0xPpY!`SVi=Z1 z-ocZRY#KySk`$gB-e%n~8LgU@W!f9GyGfLBVQz6y2&YcHO-}13V-}~bv0-|iVNusI z)|Zmtq)q8nI9jyg*i{oe-~e}<^TzYihW_EhOLWWG;TAQ7)3jNyH>WDOH7X2ygUgKO zZO4$p%Sd1<%U%W_H@6xIxWX_zG}vZ|C606} zk$3j{xEQ|+kz7nq$z;_u$wih)nq6Ul4H~WY90ha!apb(RT)^J%di=n(OA;J&QkrQLem-Vrj?@ ze*YD(1VQdxY#RNg9{x6$h_qiq`elSwgs&hBAPgc5Aq)c)RRL#T5pr*DAcY$&u=^DC zFmm5)h*y9PV@9LFZ6wzsle91ySK@S92k@T*8p2e$@I zG86oI@ItC*X2CxPN6CAEK6GXJTVPgHL(PMo6LQ5cqwYh3^`Yf<^Ldcn=Fl*0Sq$s4 z!j8K)ZI}%X1zQu>r$ub zQK;ujD=Wn%+I#6xDBW$&q_a(^LkP8H2-C1DoWk9NrD_=s7s^iX3Ph!*eWjZW<8H(B zzz%i;#z2Q!=bp!fbY|Mr=oIXt+ptWp({?%q9{gnJMjJLo5jsK)phyJPZRdiLTQdu6 zn1B+h@p8qg{ux*o%@tiGI>|j22GHVZCrh4}h*(4YHxargALz~leuCeOCi7nh)>K+5 zG<$NOTvuMd8o0xAnjhIC!Ow=j-dhxiaTa$6_eg_^!yJDcKd@ zhVSXA^uu*)+Ax*wE$@B;RxVb664jt~b@-(TPA*N7L~JU`FcvTeeVBv|Y&5+xGlfV6 zzbw2(D!~gQbPbilM3AfZw75Cf*uoV&4*M@KcX7=X*)oIx3XXiD?l5j0MflYSZ`V*B z7whm57D5BzdjMgKL$mf=VI!kij=_qWa8WH2lDSCho+3lj9Gdp#@^6s3*ozdR!Ndn) z>Nz5|fG_5?o;D=irh?bT{ow4_n{6EN2s-?s>=(tDYzv)CE1(MFcqJJPF57$anpj?$h4~{B0HyduW;l;CLdfFPmM>Q3Aj11RuV}_9`{cLoFJkFdH(4x z<|wuSx?K$gPoye#&05>LJy`U zhdTv2jy~jNc^0TJZZks&V^7%Q6y8~{?e*G8Z0E;`<2ai(e{MtQ9|g4qQQAUFtExatQKZV`GOPLa{mq*< z@4cD(FWo%SYPGZ&3XI=>=$x4tyw#dg|0cJ_CTlFle2p2OQhOAfG`Ofzord1iCfb=k z(;;UZks6B{4buo{I=9iZ&{oiNZD`Oltdc(47JW2IW1vqkjm3HtKz*jx(ID{uzd%?( z${H@3Y3VY)=%?+p16+2qDD7M%w2LNy>!IDWhbCz+O)Y8@$>w!i=^BelSA7s&n)buu zUTGqwEJjsof@Xm308Ofqr}~PWkfT951oFNo&RbuJN0IrY!;rBNItq5_ZR~o{uwRCN zh8emG3{9%Op=5S1A5cG0m;cZn6A}_DfdMig!ohjcE8Yzbk)7hRU^;*$&H2Q=c=z(J z!8Tv275KZv^HaA{VrO|6*o;s9_=so=e-Iu72_Hfj7E9qT!UIT3*dZ=OlH@6|5IJE! z57L|uV11yYAk1ShogWqPQ)H}GRjNL!QSE(|5f*EBca zDz}4Gpl^_Ce#z+#E+Dy3g9}Q|sA@F0q=K!!!Rt#k6ojMzT+z^M-JJoOgsc8)pc<^| z)lk)#hRtR)>br?KHLiYXm@p2S$IhE}!DR)DyqPG7;`<}~!!bt%rdaD{CQT348O7M#U9alA=v{jBJhMw7U zYKZCB$1Spz2pnW~+)R8YyAVTE#FFz$} z@u5H({Ta)j$3wo_38Y^}cm?4kLL8wTp#z~4Afq~f2Q>K>p=KK7K&xbq!YIao=>uR#au@I~Ns0YHQBX(b`sFb8_iD3;w)+09a$ zb3^=+{D@eh&}&{qCBK3|5w_bm>R9ub!i{b3gQbHhS7hV!w}Hy2UTCgJS-!EuZ@rAo zpghPg!;jA*6cI`QUVtswj_dF_WVC$Iq8W{IRLg*5EYfvtPng# zB9|dD6p@9kl}a>x$Uw!Vbh!5+7&=h^C1j9#tyHL3itws2hH#W^IF!oDe>Z@#Ve@TlCM7TS~{U;rNJXbeAC|@yN){yc)@Iz7IWDwpA`@L zN5Z>tF;;_jiKN*t&Y6iCUag6!Gn;1?FN~M$63g2K#cXD<*iaj-wNHl1{<*zDNu j@cGPmazCVXHImjKxsk__)_~N;N%~P_PpdBW4vhW>ef1%B diff --git a/try2/shared/classification.py b/try2/shared/classification.py index 1b80fd1..4ccb1a1 100644 --- a/try2/shared/classification.py +++ b/try2/shared/classification.py @@ -20,6 +20,7 @@ import joblib import openpyxl from sklearn.linear_model import LogisticRegression +from sklearn.ensemble import RandomForestClassifier from sklearn.preprocessing import StandardScaler from sklearn.pipeline import Pipeline @@ -96,24 +97,43 @@ def extract_ethmoid_features( volume: np.ndarray, mask: np.ndarray, class_id: int ) -> np.ndarray: """ - 6 CT intensity features sampled from voxels belonging to class_id. - volume shape: (H, W, D), dtype typically int16 (HU values). + Ethmoid filling features derived from the middle sagittal slab. + + Best config from sweep: slab ±1 voxel around the middle of the W range + (not the centroid), threshold -200 HU for filled_ratio. + + volume shape: (H, W, D) — W = columns = left-right axis. + Sagittal plane = H×D at fixed W. Returns zeros if class_id is absent. """ - voxels = np.argwhere(mask == class_id) + voxels = np.argwhere(mask == class_id) # (N, 3): h, w, d if len(voxels) == 0: return np.zeros(6, dtype=np.float32) - intensities = volume[voxels[:, 0], voxels[:, 1], voxels[:, 2]].astype(np.float32) + # Middle of the W extent of the mask (not centroid — midpoint of range) + w_min, w_max = int(voxels[:, 1].min()), int(voxels[:, 1].max()) + w_mid = (w_min + w_max) // 2 + w_lo = max(0, w_mid - 1) + w_hi = min(volume.shape[1] - 1, w_mid + 1) + + slab_mask = (voxels[:, 1] >= w_lo) & (voxels[:, 1] <= w_hi) + slab_voxels = voxels[slab_mask] + + if len(slab_voxels) == 0: + return np.zeros(6, dtype=np.float32) + + intensities = volume[ + slab_voxels[:, 0], slab_voxels[:, 1], slab_voxels[:, 2] + ].astype(np.float32) - mean_hu = float(np.mean(intensities)) - median_hu = float(np.median(intensities)) - std_hu = float(np.std(intensities)) - p10_hu = float(np.percentile(intensities, 10)) - p90_hu = float(np.percentile(intensities, 90)) - # fraction of voxels above air threshold — fluid/mucus >> -500 HU - filled_ratio = float(np.mean(intensities > -500)) + mean_hu = float(np.mean(intensities)) + median_hu = float(np.median(intensities)) + std_hu = float(np.std(intensities)) + p10_hu = float(np.percentile(intensities, 10)) + p90_hu = float(np.percentile(intensities, 90)) + # -200 HU threshold: better air/soft-tissue boundary than -500 + filled_ratio = float(np.mean(intensities > -200)) return np.array( [mean_hu, median_hu, std_hu, p10_hu, p90_hu, filled_ratio], @@ -182,7 +202,14 @@ def build_feature_matrix( return X_dict, y_dict -def _make_clf() -> Pipeline: +def _make_clf(task: str) -> Pipeline: + """Return the best classifier pipeline for the given task. + + All tasks use LogisticRegression — stable under LOO-CV with small + positive class sizes (23–26 positives out of 133). RandomForest showed + higher variance and collapsed recall in LOO-CV despite appearing better + in a fixed-split sweep. + """ return Pipeline([ ("scaler", StandardScaler()), ("clf", LogisticRegression(class_weight="balanced", max_iter=1000, C=1.0)), @@ -208,7 +235,7 @@ def train_classifiers( y = y_dict[task] if len(X) == 0: continue - clf = _make_clf() + clf = _make_clf(task) clf.fit(X, y) clfs[task] = clf return clfs From 43aa36994537b44a1732b224ad78bcf8cc77aca9 Mon Sep 17 00:00:00 2001 From: Petra Hedesiu Date: Thu, 28 May 2026 16:30:33 -0400 Subject: [PATCH 5/5] merged left and right sides --- try2/experiments/exp_classification.py | 70 ++++++---- .../classification/ext_predictions.csv | 59 +++++--- .../__pycache__/classification.cpython-39.pyc | Bin 8537 -> 9048 bytes try2/shared/classification.py | 131 ++++++++++-------- 4 files changed, 155 insertions(+), 105 deletions(-) diff --git a/try2/experiments/exp_classification.py b/try2/experiments/exp_classification.py index 69d2ce4..a9e3f8d 100644 --- a/try2/experiments/exp_classification.py +++ b/try2/experiments/exp_classification.py @@ -46,6 +46,7 @@ predict_classifications, TASK_TO_FEATURES, _make_clf, + extract_features, ) from shared.config import ExperimentConfig @@ -159,20 +160,27 @@ def load_ext_data(ext_dir: str): # LOO cross-validation # --------------------------------------------------------------------------- -def loo_evaluate(X: np.ndarray, y: np.ndarray, task: str) -> dict: - if len(X) < 4: - print(f" [skip] {task}: only {len(X)} samples with matched labels") +def loo_evaluate( + X: np.ndarray, y: np.ndarray, patient_idx: np.ndarray, task: str +) -> dict: + """Patient-level LOO-CV: hold out both sides of one patient per fold.""" + patients = np.unique(patient_idx) + if len(patients) < 4: + print(f" [skip] {task}: only {len(patients)} patients with matched labels") return {} - loo = LeaveOneOut() y_true, y_pred, y_prob = [], [], [] - for train_idx, test_idx in loo.split(X): + for p in patients: + test_mask = patient_idx == p + train_mask = ~test_mask clf = _make_clf(task) - clf.fit(X[train_idx], y[train_idx]) - y_true.append(y[test_idx[0]]) - y_pred.append(clf.predict(X[test_idx])[0]) - y_prob.append(clf.predict_proba(X[test_idx])[0, 1]) + clf.fit(X[train_mask], y[train_mask]) + probs = clf.predict_proba(X[test_mask])[:, 1] + preds = clf.predict(X[test_mask]) + y_true.extend(y[test_mask].tolist()) + y_pred.extend(preds.tolist()) + y_prob.extend(probs.tolist()) y_true = np.array(y_true) y_pred = np.array(y_pred) @@ -209,23 +217,28 @@ def run_ext_validation( print("No EXT patients loaded — check ext_dir structure.") return - # Run predictions + # Run predictions — each patient produces left + right rows + side_labels = { + "left": ("aeal_roof_contact", "left_ethmoid_filled"), + "right": ("aear_roof_contact", "right_ethmoid_filled"), + } rows = [] for code, vol, msk in zip(codes, volumes, masks): preds = predict_classifications(vol, msk, clfs) - row = {"patient_code": code} - row.update(preds) - # Attach ground-truth labels if available in Excel - if code in labels: - for task in TASK_TO_FEATURES: - row[f"{task}_gt"] = labels[code][task] - rows.append(row) + for side, (roof_lbl, eth_lbl) in side_labels.items(): + row = {"patient_code": code, "side": side} + row["roof_contact"] = preds[side].get("roof_contact", "") + row["ethmoid_filled"] = preds[side].get("ethmoid_filled", "") + if code in labels: + row["roof_contact_gt"] = labels[code][roof_lbl] + row["ethmoid_filled_gt"] = labels[code][eth_lbl] + rows.append(row) # Save predictions CSV os.makedirs(output_dir, exist_ok=True) csv_path = os.path.join(output_dir, "ext_predictions.csv") - fieldnames = ["patient_code"] + list(TASK_TO_FEATURES.keys()) + \ - [f"{t}_gt" for t in TASK_TO_FEATURES] + fieldnames = ["patient_code", "side", "roof_contact", "ethmoid_filled", + "roof_contact_gt", "ethmoid_filled_gt"] with open(csv_path, "w", newline="") as f: writer = csv.DictWriter(f, fieldnames=fieldnames, extrasaction="ignore") writer.writeheader() @@ -233,13 +246,13 @@ def run_ext_validation( print(f"Predictions saved → {csv_path}") # Compute metrics for patients that have Excel labels - labeled_rows = [r for r in rows if f"{list(TASK_TO_FEATURES)[0]}_gt" in r] + labeled_rows = [r for r in rows if "roof_contact_gt" in r] if not labeled_rows: print("No EXT patients found in Excel — predictions saved, no metrics computed.") return - print(f"\nMetrics on {len(labeled_rows)} EXT patients with Excel labels:") - for task in TASK_TO_FEATURES: + print(f"\nMetrics on {len(labeled_rows)} EXT instances with Excel labels:") + for task in ("roof_contact", "ethmoid_filled"): y_true = np.array([r[f"{task}_gt"] for r in labeled_rows]) y_pred = np.array([r[task] for r in labeled_rows]) if len(np.unique(y_true)) < 2: @@ -293,13 +306,14 @@ def main(): print("\nERROR: No patients matched. Check patient code extraction.") return - X_dict, y_dict = build_feature_matrix(codes, volumes, masks, labels) - n = len(y_dict["aeal_roof_contact"]) - print(f"\nFeature matrix built: {n} patients with matched labels") + X_dict, y_dict, patient_idx = build_feature_matrix(codes, volumes, masks, labels) + n_patients = len(np.unique(patient_idx)) + n_samples = len(y_dict["roof_contact"]) + print(f"\nFeature matrix built: {n_patients} patients → {n_samples} samples (L+R pooled)") - # LOO cross-validation + # LOO cross-validation (patient-level) print(f"\n{'=' * 60}") - print("Leave-One-Out Cross-Validation (CROP1)") + print("Patient-level Leave-One-Out Cross-Validation (CROP1)") print(f"{'=' * 60}") for task, feat_key in TASK_TO_FEATURES.items(): @@ -308,7 +322,7 @@ def main(): print(f"\n--- {task} ---") print(f" Features : {feat_key} shape={X.shape}") print(f" Positive : {int(y.sum())}/{len(y)} ({100*y.mean():.1f}%)") - res = loo_evaluate(X, y, task) + res = loo_evaluate(X, y, patient_idx, task) if res: print(f" Accuracy : {res['accuracy']:.3f}") print(f" Balanced accuracy : {res['balanced_accuracy']:.3f}") diff --git a/try2/outputs/classification/ext_predictions.csv b/try2/outputs/classification/ext_predictions.csv index 92cc202..321b3d5 100644 --- a/try2/outputs/classification/ext_predictions.csv +++ b/try2/outputs/classification/ext_predictions.csv @@ -1,20 +1,39 @@ -patient_code,aeal_roof_contact,aear_roof_contact,left_ethmoid_filled,right_ethmoid_filled,aeal_roof_contact_gt,aear_roof_contact_gt,left_ethmoid_filled_gt,right_ethmoid_filled_gt -EXT01,1,0,1,0,,,, -EXT03,0,0,0,1,,,, -EXT04,0,1,1,1,,,, -EXT05,0,0,0,1,,,, -EXT06,1,1,1,1,,,, -EXT07,0,1,1,0,,,, -EXT08,0,1,0,0,,,, -EXT09,0,0,0,0,,,, -EXT10,0,0,1,1,,,, -EXT11,0,1,1,1,,,, -EXT12,0,0,0,1,,,, -EXT13,1,0,1,1,,,, -EXT14,0,1,0,1,,,, -EXT15,1,1,1,0,,,, -EXT17,0,0,0,1,,,, -EXT18,1,1,1,1,,,, -EXT19,1,1,1,0,,,, -EXT24,1,1,0,1,,,, -EXT25,0,0,1,1,,,, +patient_code,side,roof_contact,ethmoid_filled,roof_contact_gt,ethmoid_filled_gt +EXT01,left,1,1,, +EXT01,right,0,1,, +EXT03,left,0,0,, +EXT03,right,0,1,, +EXT04,left,0,1,, +EXT04,right,0,1,, +EXT05,left,0,0,, +EXT05,right,0,1,, +EXT06,left,1,1,, +EXT06,right,0,1,, +EXT07,left,0,1,, +EXT07,right,0,1,, +EXT08,left,0,0,, +EXT08,right,1,0,, +EXT09,left,0,0,, +EXT09,right,0,0,, +EXT10,left,0,1,, +EXT10,right,0,1,, +EXT11,left,0,1,, +EXT11,right,1,1,, +EXT12,left,0,0,, +EXT12,right,0,1,, +EXT13,left,1,1,, +EXT13,right,0,1,, +EXT14,left,0,0,, +EXT14,right,0,1,, +EXT15,left,1,1,, +EXT15,right,1,0,, +EXT17,left,0,0,, +EXT17,right,0,1,, +EXT18,left,1,1,, +EXT18,right,1,1,, +EXT19,left,1,1,, +EXT19,right,1,0,, +EXT24,left,1,1,, +EXT24,right,1,1,, +EXT25,left,0,1,, +EXT25,right,0,1,, diff --git a/try2/shared/__pycache__/classification.cpython-39.pyc b/try2/shared/__pycache__/classification.cpython-39.pyc index 25935418b77cc9c5364f33a653652286126ad7b2..47347a9a68ba4b9f48a38499ef86a0daabe10abf 100644 GIT binary patch delta 3830 zcmcIm-)|g89pBmA+dJQ#owMV_j$^xHJC1#clhzIda1zwgxU@7)N+VD(Nl)vWIp5mr z-Sh08e>8g>D3Pj4AW^owStv$AA`${sdEkv_UV#@_A&@{+$g1*YJ%%56s7Ixck?Ho?WkCqL3}k(Aj@6eaw7YOHw? zBSEQLJFFW*ST$w*E?~rxcuqsqY@2$%D`2@EubIO1xxGz2mu|v}O-F7;x`Jue;vIFw z{2&IFLzn874dLacf2A&XY0ZRCZNZ`7N+>nCX_wj$X9kBv=Md>CP^T?DYcxWB(&EQs?cnnAdT@FK!X z2y+OR5iTKIK$u1N9Kz?rFBC=w=Al{a6KyL(*_Cj$@UUZ9|1OLW77q4a>wXn2w3bD< ze(FTH+518EESll^iea|z+?&_1D>cuM896(tINIt#LFyMz1dDr1yEyGt;8hRtuY8j36veL--icbP1Ac6WodqA!$ z_kMX-ZFJq$8rep!kq@$gwn4

(jm(kTF~8O?CGsu26|mmOS)jG4i83$hjfRB%ux^P&FMfFpZNNtRG{w?Cj%{m zf%72sLq+V#@l^Xb%kHR+W4^wU4h+By%l$x=y>!xJdC({X$AW?!2amBzxv_B!Ob?{_VRT(fi|=yyNq9dhX^7y$m!r-XhS$v4-?S z>Gt!D6IyN=#4eictei*^N3D&OEVdP|0w*tGc=$4qs40ZG4H2o;fMY2~riM);gIoRq zrGhr&Tk$4wGaa9Acukd!BB;nk8nyqJ+q zm~_$Djz?14CFn-Ewv$*))$l+B4S#m>oOy(J`VU`$!$)8j=JAj+l|q(eiLQ>43{jz# zQF}=I%^tM!+DS6FXQW_!&(OuYPc9GqgWM>Aam?d#!g1(EUGT=N@>1EZ{2sbq20+QE z;JUCUD$o`;5Z*$#2>{G`7u)xMS^qpP2_3Z}tz!AbOiP9XI(YO%DRhTUPMK(NWRJZ- zt(F7xCNH}8N%)JQk<1$C!XzueMOp9gs6YB_dgbjem;9HjGd&0wxZo4Kc>S7DUa(X$l)IPG00;;yp z{4D&{@QLPqumlH(@0Ro;UkD9E$)^h%7;^bS`~(`lg6{w@yv7Jri5jq=ZYd4m)kkDk zftO(mUIqhx=|;u^x(|%23ey91LlG5t9aIMN4#fU#Edb7bq_Y$=cGX?Np!WdiH~?N} zI+$dn36&h;$^(LmiLi~!UV3r%kPKVpCFsF$W2Qxj3FnMMiP)_Ezg}WCe5@Wa4YUs! z21rB4r9~-Oe~Mdh6tbpy3N_4r%sFqvSTQBZLR71$fY}(OTlrPLreffX@|)tRvgrv@fCU(H6+7IIzrwPn+Ql9@-X8%7$Y-CRl12Xyb z(LIfwWAuG}d4y}chj*9NGYaovDjNe!fsKP@6!b?nHnFF&N%lk-k)f_$^k!3gI-6!i z@TB2i$iSJgBkhyyJZK9{eO)QeY%Yh_)t|pUM^~JRE$Y;@SGeoZ8W*%?`wp-8^O+2# zOVw(bmZ}xsF8OqNl~*^oFPswn6D7CpQECI2sR;<{*1oo$Iz&H#3+D6r3)N?AE9`$_NhrelIGpFiSSDljW z${}C;E_EtCufT|PXvI?yd=-YrksHQ99!GY0*{3`X8U8>U*C{+9P50)(D4f-`6V`{X z>9@!?cAsD3{-&r<=Di_c|MCWu>B?2ddKDBIj|PH#iY^!uJKf{)79{^{QD>7`=hQAawq&v z*YoW!VHXX$V2|)$T_0tpu!)m?5dO0J$ycwTD81Cg3l+x}4`P((JHCM9W|cXd0>Qlx zLF>s1GIB(|hvJp6qi3>t4!dUoN;*L6>48QXelNn$^NFIY!mBC`y^%^olQp0CEn8q+0z?czIN9Lpx50r%l=rE zQrihB6HAkcG21e-9CYsvx&k0BAer_}7qTDD_6;?WFWxOk23Rc2)9G6?bHy1NE0&Ha zlA^T|rz_~x(j|UZiVYwPeCn30)i^d#jC+D&;;0zF=EYO=)3POj;&6zC@I9j^(=vky zq|r>oc-Yr~a>?CfoIOW{Bp8qSYupiX z=(6MbTmaKJ6<7-Qs0|^x>t2vnaUD*uc)5vq1L0c;c)Ud~LLWjuKv5NtXhn1)EFd6f zh;D$Q8ELytO}ql_$iO*Pm?VrS1uR=yk5YW2=GSGpM=2oJ%a@8h;#KUvh;UO;nu@?e zTh3xZ#wlJrE-Wo)qy_fLT0GWXTF%z3qisoq_%pDasar_W7W6@bpcsgbYzHp)qY|U71|xiYbO7Bt)KLh zK{9w~q!ZhrVTiXMemd}X@&=OYXw0{Zm0QH2Ozc+3|~JVnc-vu|2UIo*U|)x{WOxM~CSG$LO@qABIP|)3~&5FA4&@6|SBTMg7>&T<)@sU3Zt&zUw7W9Mw z@z>)hiez^U#MXm}Pr*;J`&-b6SsPecbAaQ0vacQ}z~5JItBq8{2vUIx>^*(2KhOdo z*nN!=rXHw`w66zh00g|B3V;JyAMo{k66itdssheB@Rz>(L3nGl*sMW?!leMp6t6rp zNgdNx?n-c9s#avhCY~V5QnusE1YDus^`WO#u{ZcUy=Ya1HQPvP9KY+2;5+_UmT=LK zl^W-}1g0o07P+@(!WOQ{Ridxwq-UrpWAT{Ij!ap6zL3g`Qe7HphaFp_I=JW5*@tClQFxZyv>rm`~_1kbe?2uBJ|a|oN5 ze=z|Nr7O%9!mf*>3r)OhZ=;7K16eOvm_i`n*_OGDZ(fqoznr6Nje|>MUfrm&O_yJT z`6?ckt=?&p!8cVCw#bYxKP0Eo?z5L+9be>Tb=bb#=eZ6;+Plf Dict[str, Dict[str, int]]: - """Return dict keyed by patient_code with 4 binary label values.""" + """Return dict keyed by patient_code. + + Each entry has per-side labels (still stored separately so we can pool + them into side-agnostic samples in build_feature_matrix): + aeal_roof_contact, aear_roof_contact, + left_ethmoid_filled, right_ethmoid_filled + """ wb = openpyxl.load_workbook(xlsx_path) ws = wb.active labels: Dict[str, Dict[str, int]] = {} @@ -159,47 +163,48 @@ def build_feature_matrix( volumes: List[np.ndarray], masks: List[np.ndarray], labels: Dict[str, Dict[str, int]], -) -> Tuple[Dict[str, np.ndarray], Dict[str, np.ndarray]]: +) -> Tuple[Dict[str, np.ndarray], Dict[str, np.ndarray], np.ndarray]: """ + Pool left and right sides into side-agnostic feature matrices. + + Each patient contributes 2 rows (left side first, then right side). Returns: - X_dict — {"roof_L": (N,4), "roof_R": (N,4), "ethmoid_L": (N,6), "ethmoid_R": (N,6)} - y_dict — {"aeal_roof_contact": (N,), "aear_roof_contact": (N,), - "left_ethmoid_filled": (N,), "right_ethmoid_filled": (N,)} - Only patients whose code appears in `labels` are included; returns their - filtered indices as well. + X_dict — {"roof": (2N, 6), "ethmoid": (2N, 6)} + y_dict — {"roof_contact": (2N,), "ethmoid_filled": (2N,)} + patient_idx — (2N,) integer array mapping each row to its patient index. + Used for patient-level LOO-CV (hold out both sides together). """ - feat_roof_L, feat_roof_R = [], [] - feat_eth_L, feat_eth_R = [], [] - y_aeal_roof, y_aear_roof = [], [] - y_left_eth, y_right_eth = [], [] + feat_roof, feat_eth = [], [] + y_roof, y_eth = [], [] + patient_idx = [] - for code, vol, msk in zip(patient_codes, volumes, masks): + for p_idx, (code, vol, msk) in enumerate(zip(patient_codes, volumes, masks)): if code not in labels: continue feats = extract_features(vol, msk) - feat_roof_L.append(feats["roof_L"]) - feat_roof_R.append(feats["roof_R"]) - feat_eth_L.append(feats["ethmoid_L"]) - feat_eth_R.append(feats["ethmoid_R"]) - lbl = labels[code] - y_aeal_roof.append(lbl["aeal_roof_contact"]) - y_aear_roof.append(lbl["aear_roof_contact"]) - y_left_eth.append(lbl["left_ethmoid_filled"]) - y_right_eth.append(lbl["right_ethmoid_filled"]) + lbl = labels[code] + # Left side (class_id=1) + feat_roof.append(feats["roof_L"]) + feat_eth.append(feats["ethmoid_L"]) + y_roof.append(lbl["aeal_roof_contact"]) + y_eth.append(lbl["left_ethmoid_filled"]) + patient_idx.append(p_idx) + # Right side (class_id=2) + feat_roof.append(feats["roof_R"]) + feat_eth.append(feats["ethmoid_R"]) + y_roof.append(lbl["aear_roof_contact"]) + y_eth.append(lbl["right_ethmoid_filled"]) + patient_idx.append(p_idx) X_dict = { - "roof_L": np.stack(feat_roof_L) if feat_roof_L else np.empty((0, 6)), - "roof_R": np.stack(feat_roof_R) if feat_roof_R else np.empty((0, 6)), - "ethmoid_L": np.stack(feat_eth_L) if feat_eth_L else np.empty((0, 6)), - "ethmoid_R": np.stack(feat_eth_R) if feat_eth_R else np.empty((0, 6)), + "roof": np.stack(feat_roof) if feat_roof else np.empty((0, 6)), + "ethmoid": np.stack(feat_eth) if feat_eth else np.empty((0, 6)), } y_dict = { - "aeal_roof_contact": np.array(y_aeal_roof, dtype=np.int32), - "aear_roof_contact": np.array(y_aear_roof, dtype=np.int32), - "left_ethmoid_filled": np.array(y_left_eth, dtype=np.int32), - "right_ethmoid_filled": np.array(y_right_eth, dtype=np.int32), + "roof_contact": np.array(y_roof, dtype=np.int32), + "ethmoid_filled": np.array(y_eth, dtype=np.int32), } - return X_dict, y_dict + return X_dict, y_dict, np.array(patient_idx, dtype=np.int32) def _make_clf(task: str) -> Pipeline: @@ -217,10 +222,8 @@ def _make_clf(task: str) -> Pipeline: TASK_TO_FEATURES = { - "aeal_roof_contact": "roof_L", - "aear_roof_contact": "roof_R", - "left_ethmoid_filled": "ethmoid_L", - "right_ethmoid_filled": "ethmoid_R", + "roof_contact": "roof", + "ethmoid_filled": "ethmoid", } @@ -260,18 +263,26 @@ def predict_classifications( volume: np.ndarray, mask: np.ndarray, clfs: Dict[str, Pipeline], -) -> Dict[str, int]: +) -> Dict[str, Dict[str, int]]: """ - Given a single patient's volume and predicted mask, return classification - predictions as a dict of {task_name: 0_or_1}. + Given a single patient's volume and mask, return side-agnostic predictions. + + Returns: {"left": {"roof_contact": 0/1, "ethmoid_filled": 0/1}, + "right": {"roof_contact": 0/1, "ethmoid_filled": 0/1}} """ feats = extract_features(volume, mask) - results: Dict[str, int] = {} - for task, feat_key in TASK_TO_FEATURES.items(): - if task not in clfs: - continue - x = feats[feat_key].reshape(1, -1) - results[task] = int(clfs[task].predict(x)[0]) + side_feat = {"left": ("roof_L", "ethmoid_L"), "right": ("roof_R", "ethmoid_R")} + results: Dict[str, Dict[str, int]] = {} + for side, (rf_key, eth_key) in side_feat.items(): + results[side] = {} + if "roof_contact" in clfs: + results[side]["roof_contact"] = int( + clfs["roof_contact"].predict(feats[rf_key].reshape(1, -1))[0] + ) + if "ethmoid_filled" in clfs: + results[side]["ethmoid_filled"] = int( + clfs["ethmoid_filled"].predict(feats[eth_key].reshape(1, -1))[0] + ) return results @@ -279,13 +290,19 @@ def predict_classifications_proba( volume: np.ndarray, mask: np.ndarray, clfs: Dict[str, Pipeline], -) -> Dict[str, float]: +) -> Dict[str, Dict[str, float]]: """Same as predict_classifications but returns probability of positive class.""" feats = extract_features(volume, mask) - results: Dict[str, float] = {} - for task, feat_key in TASK_TO_FEATURES.items(): - if task not in clfs: - continue - x = feats[feat_key].reshape(1, -1) - results[task] = float(clfs[task].predict_proba(x)[0, 1]) + side_feat = {"left": ("roof_L", "ethmoid_L"), "right": ("roof_R", "ethmoid_R")} + results: Dict[str, Dict[str, float]] = {} + for side, (rf_key, eth_key) in side_feat.items(): + results[side] = {} + if "roof_contact" in clfs: + results[side]["roof_contact"] = float( + clfs["roof_contact"].predict_proba(feats[rf_key].reshape(1, -1))[0, 1] + ) + if "ethmoid_filled" in clfs: + results[side]["ethmoid_filled"] = float( + clfs["ethmoid_filled"].predict_proba(feats[eth_key].reshape(1, -1))[0, 1] + ) return results