{"cells":[{"cell_type":"markdown","metadata":{},"source":"# Role 2 preprocessing v2 — shard 15/23\n\nuint16 raw volumes, correct slice ordering (projection onto slice normal),\nper-plane laterality canonicalization, failures.csv audit trail."},{"cell_type":"code","metadata":{},"source":"%%writefile role2_pipeline.py\n\"\"\"\nRole 2 preprocessing pipeline v2 (fixed).\n\nFixes vs v1, each tied to evidence:\n\n1. Slice ordering. v1 sorted by ImagePositionPatient[2] (z) alone. The\n   full-dataset scan showed ~18% of series (2,965 sagittal + 1,431 coronal,\n   ~139k slices) have |slice-normal z-component| ~ 0, so every sort key\n   ties and Python's stable sort preserves arbitrary os.listdir() order.\n   Fix: sort by position projected onto the slice normal,\n   key = dot(ImagePositionPatient, cross(row_dir, col_dir)), which is\n   correct for ANY gantry tilt; falls back to InstanceNumber.\n\n2. Laterality. v1 always did np.flip(volume, axis=-1). That is a\n   left-right mirror only for axial/coronal. For sagittal series the\n   left-right axis IS the through-plane axis, so the anatomically correct\n   mirror is reversing slice order -- v1 was flipping anterior-posterior\n   instead, producing impossible anatomy. Fix: derive the left-right\n   carrying axis from ImageOrientationPatient per plane.\n   (Flip is still gated on the Laterality tag, per ROLES.md's\n   \"laterality_flipped must be auditable\".)\n\n3. Storage. v1 stored float32 with a baked-in 0.5/99.5 percentile\n   normalization: ~75 MB/study, irreversible. Fix: uint16 raw after\n   slope/intercept rescale, clipped to [0, 65535] (~35 MB/study) --\n   normalization belongs to the training DataLoader where it is\n   reversible and tunable per experiment.\n\n4. Silent drops. v1 caught per-series exceptions, printed [WARN], and\n   still marked the study done. Fix: every failed series is recorded in\n   <out>/failures.csv, and a study is only marked done when every series\n   train_series.csv expects has a disposition (written or failed).\n\n5. Shards instead of one dataset per batch. Deterministic assignment\n   sha1(StudyInstanceUID) % n_shards, so a study always lands in the same\n   shard no matter which session processes it. Shards are created once\n   and mounted by Role 3 (ROLES.md contract: no ad hoc upload/download\n   between stages).\n\"\"\"\n\nimport hashlib\nimport os\nfrom datetime import datetime, timezone\n\nimport numpy as np\nimport pandas as pd\nimport pydicom\nimport h5py\n\n# ---------------------------------------------------------------- geometry\n\ndef slice_normal(ds):\n    \"\"\"Through-plane direction from ImageOrientationPatient, or None.\"\"\"\n    if not hasattr(ds, \"ImageOrientationPatient\"):\n        return None\n    rd = np.asarray(ds.ImageOrientationPatient[:3], dtype=np.float64)\n    cd = np.asarray(ds.ImageOrientationPatient[3:], dtype=np.float64)\n    n = np.cross(rd, cd)\n    norm = np.linalg.norm(n)\n    return n / norm if norm > 1e-9 else None\n\n\ndef canonical_normal(normal):\n    \"\"\"Sign-canonicalized slice normal.\n\n    dot(IPP, normal) ordering is correct within a series, but the normal's\n    SIGN is scanner-dependent: two scanners imaging the same knee can\n    traverse it in opposite directions. Flipping the sign so the dominant\n    component is positive makes traversal anatomically canonical and\n    scanner-independent: sagittal medial->lateral, axial feet->head,\n    coronal anterior->posterior. Any order-aware model downstream then\n    sees consistent geometry across all sites.\n    \"\"\"\n    if normal is None:\n        return None\n    if normal[np.argmax(np.abs(normal))] < 0:\n        return -normal\n    return normal\n\n\ndef sort_key(ds, normal):\n    \"\"\"Position projected onto the slice normal. Falls back to\n    InstanceNumber, then 0 (file order), when geometry is unavailable.\"\"\"\n    if normal is not None and hasattr(ds, \"ImagePositionPatient\"):\n        return float(np.dot(np.asarray(ds.ImagePositionPatient, dtype=np.float64), normal))\n    return None\n\n\n# ------------------------------------------------------------- laterality\n\nCANONICAL_LATERALITY = \"L\"\n# DICOM patient coords: +x = patient Left, +y = posterior, +z = head.\n\ndef planar_flip_info(normal, first_ds, plane_csv=None):\n    \"\"\"Decide which array axis carries the left-right direction and how to\n    mirror a series onto the canonical (L) chirality.\n\n    Priority:\n      1. Physical ImageOrientationPatient normal vector (ground truth).\n      2. If geometry is unavailable, fallback to train_series.csv plane label.\n      3. Return (\"none\", None) if orientation is ambiguous.\n\n    Returns (method, axis_or_none):\n      (\"slice-reverse\", None)  -- sagittal-type: L-R is through-plane,\n                                  canonicalize by reversing slice order\n      (\"in-plane\", axis)       -- axial/coronal-type: flip that image axis\n      (\"none\", None)           -- no reliable L-R axis derivable\n    \"\"\"\n    # 1. Primary: Ground truth from ImageOrientationPatient normal vector\n    if normal is not None:\n        dominant_axis = int(np.argmax(np.abs(normal)))\n        if dominant_axis == 0:\n            # Dominant component is Left-Right (X) -> Sagittal-type (includes oblique up to ~55 deg)\n            return (\"slice-reverse\", None)\n\n        # Dominant normal is Y (Anterior-Posterior / Coronal) or Z (Superior-Inferior / Axial)\n        if hasattr(first_ds, \"ImageOrientationPatient\"):\n            rd = np.asarray(first_ds.ImageOrientationPatient[:3], dtype=np.float64)\n            cd = np.asarray(first_ds.ImageOrientationPatient[3:], dtype=np.float64)\n            rx, cx = abs(rd[0]), abs(cd[0])\n\n            # Guard: If neither in-plane direction has meaningful L-R alignment\n            if rx < 0.2 and cx < 0.2:\n                return (\"none\", None)\n\n            if rx >= cx:\n                return (\"in-plane\", 2)  # rd aligns with image columns (numpy axis 2)\n            else:\n                return (\"in-plane\", 1)  # cd aligns with image rows    (numpy axis 1)\n\n    # 2. Fallback: Metadata plane hint when ImageOrientationPatient geometry is unavailable\n    if plane_csv is not None:\n        plane_clean = str(plane_csv).strip().lower()\n        if plane_clean == \"sagittal\":\n            return (\"slice-reverse\", None)\n        elif plane_clean in (\"coronal\", \"axial\"):\n            return (\"in-plane\", 2)\n\n    return (\"none\", None)\n\n\ndef apply_canonical_laterality(volume, normal, first_ds, laterality, plane_csv=None):\n    \"\"\"Mirror onto canonical 'L' chirality when the study side is 'R'.\n    Returns (volume, flipped_bool, method_str).\"\"\"\n    if laterality is None:\n        return volume, False, \"unknown-laterality\"\n    if laterality == CANONICAL_LATERALITY:\n        return volume, False, \"already-canonical\"\n    method, axis = planar_flip_info(normal, first_ds, plane_csv=plane_csv)\n    if method == \"slice-reverse\":\n        return np.ascontiguousarray(volume[::-1]), True, \"slice-reverse\"\n    if method == \"in-plane\":\n        return np.ascontiguousarray(np.flip(volume, axis=axis)), True, f\"in-plane-axis{axis}\"\n    return volume, False, \"no-lr-axis\"\n\n\n# --------------------------------------------------------------- processing\n\ndef load_series(series_dir):\n    \"\"\"Read all slices once (full read -- pixels are needed anyway) and\n    order them by through-plane position.\"\"\"\n    files = sorted(os.path.join(series_dir, f) for f in os.listdir(series_dir)\n                   if f.endswith(\".dcm\"))\n    if not files:\n        return None, \"no-dcm-files\"\n\n    datasets = [pydicom.dcmread(f) for f in files]\n    normal = canonical_normal(slice_normal(datasets[0])) if datasets else None\n    keys = [sort_key(ds, normal) for ds in datasets]\n\n    if normal is None or any(k is None for k in keys) or len(set(keys)) < len(keys):\n        # geometry missing or degenerate -> InstanceNumber, then file order\n        datasets.sort(key=lambda ds: int(getattr(ds, \"InstanceNumber\", 0)))\n        order_method = \"instance-number-fallback\"\n    else:\n        datasets = [ds for _, ds in sorted(zip(keys, datasets), key=lambda t: t[0])]\n        order_method = \"projection-on-canonical-normal\"\n    return datasets, order_method\n\n\ndef to_uint16(ds):\n    \"\"\"Rescale to true intensities and store as uint16. Negative values\n    (possible when intercept < 0) are clipped to 0 -- knee MRI signal is\n    non-negative in practice; the raw float path is gone by design.\"\"\"\n    pixels = ds.pixel_array.astype(np.float32)\n    slope = float(getattr(ds, \"RescaleSlope\", 1.0) or 1.0)\n    intercept = float(getattr(ds, \"RescaleIntercept\", 0.0) or 0.0)\n    values = pixels * slope + intercept\n    return np.clip(values, 0, 65535).astype(np.uint16)\n\n\ndef normalize_sequence_type(description: str) -> str:\n    \"\"\"Map raw SeriesDescription to canonical sequence taxonomy (T1, T2, PD, other).\"\"\"\n    desc = str(description).lower()\n    if \"pd\" in desc or \"proton\" in desc:\n        return \"PD\"\n    if \"t2\" in desc:\n        return \"T2\"\n    if \"t1\" in desc:\n        return \"T1\"\n    return \"other\"\n\n\ndef process_series(series_dir, plane_csv):\n    \"\"\"Process one series. Returns (result_dict, None) or (None, error_str).\"\"\"\n    datasets, order_method = load_series(series_dir)\n    if datasets is None:\n        return None, \"no-dcm-files\"\n\n    first = datasets[0]\n    normal = canonical_normal(slice_normal(first))\n    volume = np.stack([to_uint16(ds) for ds in datasets], axis=0)\n\n    plane_clean = str(plane_csv).strip().lower() if plane_csv is not None else \"unknown\"\n    laterality = str(getattr(first, \"Laterality\", \"\")).upper() or None\n    volume, flipped, flip_method = apply_canonical_laterality(\n        volume, normal, first, laterality, plane_csv=plane_clean)\n\n    raw_desc = str(getattr(first, \"SeriesDescription\", \"unknown\"))\n    seq_type = normalize_sequence_type(raw_desc)\n\n    result = {\n        \"array\": volume,\n        \"plane\": plane_clean,                # lowercase schema-compliant (\"sagittal\", \"coronal\", \"axial\")\n        \"sequence_type\": seq_type,          # canonical \"T1\" | \"T2\" | \"PD\" | \"other\"\n        \"sequence_type_raw\": raw_desc,      # raw SeriesDescription for audit\n        \"n_slices\": int(volume.shape[0]),\n        \"laterality_tag\": laterality or \"\",\n        \"laterality_flipped\": bool(flipped),\n        \"flip_method\": flip_method,\n        \"slice_order_method\": order_method,\n        \"array_shape\": str(tuple(int(x) for x in volume.shape)),\n        \"photometric\": str(getattr(first, \"PhotometricInterpretation\", \"\")),\n    }\n    return result, None\n\n\n# ------------------------------------------------------------------ shards\n\ndef study_shard(study_uid, n_shards):\n    \"\"\"Deterministic shard assignment, stable across sessions/machines.\"\"\"\n    return int(hashlib.sha1(study_uid.encode()).hexdigest(), 16) % n_shards\n\n\ndef _append_csv(path, rows, columns):\n    if not rows:\n        return\n    df = pd.DataFrame(rows, columns=columns)\n    df.to_csv(path, mode=\"a\", header=not os.path.exists(path), index=False)\n\n\ndef run_shard(input_root, train_csv, shard_index, n_shards, out_dir):\n    \"\"\"Process every study whose shard is `shard_index` into\n    <out_dir>/shard_{i:02d}.h5 (+ manifest and failures sidecars).\n    Resumable: studies already fully dispositioned in this shard are skipped.\n    \"\"\"\n    os.makedirs(out_dir, exist_ok=True)\n    shard_tag = f\"shard_{shard_index:02d}\"\n    h5_path = os.path.join(out_dir, f\"{shard_tag}.h5\")\n    manifest_path = os.path.join(out_dir, f\"{shard_tag}_manifest.csv\")\n    failures_path = os.path.join(out_dir, \"failures.csv\")\n    done_path = os.path.join(out_dir, f\"{shard_tag}_done_studies.txt\")\n\n    meta = pd.read_csv(train_csv)  # one row per series; authoritative metadata\n    meta[\"StudyInstanceUID\"] = meta[\"StudyInstanceUID\"].astype(str)\n    meta[\"SeriesInstanceUID\"] = meta[\"SeriesInstanceUID\"].astype(str)\n\n    all_studies = sorted(s for s in os.listdir(input_root)\n                         if os.path.isdir(os.path.join(input_root, s)))\n    my_studies = [s for s in all_studies if study_shard(s, n_shards) == shard_index]\n\n    done = set()\n    if os.path.exists(done_path):\n        with open(done_path) as f:\n            done = {line.strip() for line in f if line.strip()}\n    todo = [s for s in my_studies if s not in done]\n    print(f\"Shard {shard_index}/{n_shards}: {len(my_studies)} studies, \"\n          f\"{len(done)} already done, {len(todo)} to process\", flush=True)\n\n    manifest_cols = [\"StudyInstanceUID\", \"SeriesInstanceUID\", \"plane\",\n                     \"plane_derived\", \"fluid_sensitive\", \"fat_suppression\",\n                     \"sequence_type\", \"sequence_type_raw\", \"n_slices\", \"laterality_tag\",\n                     \"laterality_flipped\", \"flip_method\",\n                     \"slice_order_method\", \"array_shape\", \"dtype\",\n                     \"photometric\"]\n    fail_cols = [\"timestamp\", \"StudyInstanceUID\", \"SeriesInstanceUID\",\n                 \"expected_by_csv\", \"error\"]\n\n    for n, study_uid in enumerate(todo, 1):\n        study_dir = os.path.join(input_root, study_uid)\n        expected = meta[meta.StudyInstanceUID == study_uid]\n\n        with h5py.File(h5_path, \"a\") as f:  # per-study handle: a crash\n            # mid-study leaves valid datasets, not a corrupt file\n            pending_fail = []\n            pending_rows = []\n            csv_series = set(expected.SeriesInstanceUID)\n            disk_series = {s for s in os.listdir(study_dir)\n                           if os.path.isdir(os.path.join(study_dir, s))}\n            for series_uid in sorted(disk_series - csv_series):\n                # v2 iterates by csv rows (authoritative); anything on disk\n                # the csv doesn't know about is surfaced, never silently kept\n                pending_fail.append([_now(), study_uid, series_uid, 0,\n                                     \"on-disk-not-in-train_series_csv\"])\n            for _, row in expected.iterrows():\n                series_uid = row.SeriesInstanceUID\n                series_dir = os.path.join(study_dir, series_uid)\n                if f\"{study_uid}/{series_uid}\" in f:\n                    continue  # already written by an interrupted run\n                if not os.path.isdir(series_dir):\n                    pending_fail.append([_now(), study_uid, series_uid, 1,\n                                         \"series-dir-not-found\"])\n                    continue\n                try:\n                    result, err = process_series(series_dir, row.Anatomical_Plane)\n                    if result is None:\n                        pending_fail.append([_now(), study_uid, series_uid, 1, err])\n                        continue\n                    f.create_dataset(f\"{study_uid}/{series_uid}\",\n                                     data=result[\"array\"],\n                                     compression=\"gzip\", compression_opts=4)\n                    pending_rows.append([\n                        study_uid, series_uid, result[\"plane\"],\n                        derive_plane_from_result(result), row.Fluid_Sensitive,\n                        row.Fat_Suppression, result[\"sequence_type\"],\n                        result[\"sequence_type_raw\"],\n                        result[\"n_slices\"], result[\"laterality_tag\"],\n                        result[\"laterality_flipped\"], result[\"flip_method\"],\n                        result[\"slice_order_method\"], result[\"array_shape\"],\n                        str(result[\"array\"].dtype), result[\"photometric\"]])\n                except Exception as e:\n                    pending_fail.append([_now(), study_uid, series_uid, 1,\n                                         f\"{type(e).__name__}: {e}\"])\n\n        # study is done only when every csv-expected series is dispositioned\n        _append_csv(manifest_path, pending_rows, manifest_cols)\n        _append_csv(failures_path, pending_fail, fail_cols)\n        with open(done_path, \"a\") as f:\n            f.write(study_uid + \"\\n\")\n        if n % 10 == 0 or n == len(todo):\n            print(f\"  [{n}/{len(todo)}] {study_uid[:24]}... \"\n                  f\"series ok={len(pending_rows)} failed={len(pending_fail)}\",\n                  flush=True)\n\n    print(f\"Shard {shard_index} complete.\", flush=True)\n\n\ndef _now():\n    return datetime.now(timezone.utc).isoformat(timespec=\"seconds\")\n\n\ndef derive_plane_from_result(result):\n    # cheap cross-check without re-reading DICOM: v1 stored the derived\n    # plane; keep that contract by deriving from the raw sequence description\n    desc = result[\"sequence_type_raw\"].lower()\n    for plane, kws in {\"sagittal\": [\"sag\"], \"coronal\": [\"cor\"],\n                       \"axial\": [\"tra\", \"ax\", \"transverse\", \"axial\"]}.items():\n        if any(kw in desc for kw in kws):\n            return plane\n    return \"unknown\"\n","outputs":[],"execution_count":null},{"cell_type":"code","metadata":{},"source":"import os, glob\n\nROOT = next(p for p in [\n    \"/kaggle/input/competitions/rsna-knee-abnormality-detection\",\n    \"/kaggle/input/rsna-knee-abnormality-detection\",\n] if os.path.isdir(p))\nprint(\"ROOT =\", ROOT)\n\nSHARD_INDEX = 15\nN_SHARDS = 24\nOUT_DIR = \"/kaggle/working/shards\"\n\nimport role2_pipeline as rp\nrp.run_shard(\n    input_root=os.path.join(ROOT, \"train_series\"),\n    train_csv=os.path.join(ROOT, \"train_series.csv\"),\n    shard_index=SHARD_INDEX,\n    n_shards=N_SHARDS,\n    out_dir=OUT_DIR,\n)","outputs":[],"execution_count":null},{"cell_type":"code","metadata":{},"source":"import h5py, pandas as pd, os\ntag = f\"shard_{SHARD_INDEX:02d}\"\nh5_path = f\"/kaggle/working/shards/{tag}.h5\"\nwith h5py.File(h5_path, \"r\") as f:\n    studies = list(f.keys())\n    n_series = sum(len(f[s].keys()) for s in studies)\n    print(f\"studies: {len(studies)}, series: {n_series}\")\n    if studies:\n        s = studies[0]\n        k0 = list(f[s].keys())[0]\n        print(\"sample shape/dtype:\", f[s][k0].shape, f[s][k0].dtype)\nprint(f\"h5 size: {os.path.getsize(h5_path)/1e9:.2f} GB\")\nmpath = f\"/kaggle/working/shards/{tag}_manifest.csv\"\nfpath = \"/kaggle/working/shards/failures.csv\"\nprint(\"manifest rows:\", len(pd.read_csv(mpath)) if os.path.exists(mpath) else 0)\nprint(\"failures rows:\", len(pd.read_csv(fpath)) if os.path.exists(fpath) else 0)\nif os.path.exists(fpath):\n    print(pd.read_csv(fpath)[\"error\"].value_counts().head(10))","outputs":[],"execution_count":null}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.10"}},"nbformat":4,"nbformat_minor":4}