{"cells":[{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"import os\nimport glob\nimport numpy as np\nimport pandas as pd\nimport pydicom\nimport cv2\nfrom multiprocessing import Pool\nfrom tqdm import tqdm\nimport argparse\n\n# Config\ndef get_kaggle_path(filename):\n    if not os.path.exists('/kaggle/input'): return f'./{filename}'\n    for root, dirs, files in os.walk('/kaggle/input'):\n        if filename in files or filename in dirs:\n            return os.path.join(root, filename)\n            \n    # If we get here, the file was NOT found anywhere.\n    print(f\"CRITICAL ERROR: Could not find '{filename}' anywhere in /kaggle/input!\")\n    print(\"Here is exactly what is currently attached to your notebook:\")\n    for root, dirs, files in os.walk('/kaggle/input'):\n        for f in files:\n            print(f\" - {os.path.join(root, f)}\")\n    raise FileNotFoundError(f\"Missing required Kaggle dataset file: {filename}. Did you forget to click 'Add Data' -> 'Competition Data' -> 'rsna-2024-lumbar-spine-mri' in the top right of the Kaggle notebook?\")\n\nDATA_DIR = get_kaggle_path(\"train_series\") # Changed from train_images\nOUT_DIR = \"/kaggle/working/train_npy_cache\"\nDESC_CSV = get_kaggle_path(\"train_series.csv\")\nTARGET_SIZE = (256, 256)\nNUM_SLICES = 24\n\ndef process_dicom_volume(dicom_file_paths, laterality=\"R\", target_size=(256, 256), num_save_slices=24):\n    if len(dicom_file_paths) == 0:\n        return np.zeros((num_save_slices, 3, target_size[0], target_size[1]), dtype=np.uint8)\n\n    slices_data = []\n    for f in dicom_file_paths:\n        try:\n            dcm = pydicom.dcmread(f, stop_before_pixels=False)\n            pos = np.array(dcm.ImagePositionPatient, dtype=float)\n            ori = np.array(dcm.ImageOrientationPatient, dtype=float)\n            r_x, r_y = ori[:3], ori[3:]\n            normal = np.cross(r_x, r_y)\n            normal = normal / (np.linalg.norm(normal) + 1e-8)\n            depth = np.dot(pos, normal)\n            \n            pixel_array = dcm.pixel_array.astype(float) * float(getattr(dcm, 'RescaleSlope', 1.0)) + float(getattr(dcm, 'RescaleIntercept', 0.0))\n            slices_data.append((depth, pixel_array))\n        except Exception:\n            continue\n            \n    if len(slices_data) == 0:\n        return np.zeros((num_save_slices, 3, target_size[0], target_size[1]), dtype=np.uint8)\n\n    slices_data.sort(key=lambda x: x[0])\n    volume_3d = np.stack([s[1] for s in slices_data], axis=0)\n\n    if str(laterality).upper().startswith(\"L\"):\n        volume_3d = volume_3d[:, :, ::-1]\n\n    max_proj = np.max(volume_3d, axis=0)\n    thresh = max(10.0, np.percentile(max_proj, 15))\n    tissue_mask = (max_proj > thresh).astype(np.uint8)\n    \n    y_idx, x_idx = np.where(tissue_mask > 0)\n    if len(y_idx) > 0 and len(x_idx) > 0:\n        ymin, ymax = np.min(y_idx), np.max(y_idx)\n        xmin, xmax = np.min(x_idx), np.max(x_idx)\n        \n        h_box, w_box = ymax - ymin, xmax - xmin\n        ymin = max(0, ymin - int(h_box * 0.12))\n        ymax = min(volume_3d.shape[1], ymax + int(h_box * 0.12))\n        xmin = max(0, xmin - int(w_box * 0.12))\n        xmax = min(volume_3d.shape[2], xmax + int(w_box * 0.12))\n        \n        max_side = max(ymax - ymin, xmax - xmin)\n        cy, cx = (ymin + ymax) // 2, (xmin + xmax) // 2\n        half = max_side // 2\n        ymin = max(0, cy - half); ymax = min(volume_3d.shape[1], cy + half)\n        xmin = max(0, cx - half); xmax = min(volume_3d.shape[2], cx + half)\n        volume_3d = volume_3d[:, ymin:ymax, xmin:xmax]\n\n    orig_s = volume_3d.shape[0]\n    sample_indices = np.linspace(0, orig_s - 1, num_save_slices, dtype=int)\n    \n    processed_24 = []\n    for idx in sample_indices:\n        s_curr = volume_3d[idx]\n        s_next = volume_3d[min(idx + 1, orig_s - 1)]\n        \n        p50, p99_5 = np.percentile(s_curr, (50, 99.5))\n        ch0 = np.clip((s_curr - p50) / (p99_5 - p50 + 1e-5), 0, 1) * 255.0\n        \n        p5, p90 = np.percentile(s_curr, (5, 90.0))\n        ch1 = np.clip((s_curr - p5) / (p90 - p5 + 1e-5), 0, 1) * 255.0\n        \n        p10, p99 = np.percentile(s_next, (10, 99.0))\n        ch2 = np.clip((s_next - p10) / (p99 - p10 + 1e-5), 0, 1) * 255.0\n        \n        triplet = np.stack([ch0, ch1, ch2], axis=-1).astype(np.uint8)\n        resized_triplet = cv2.resize(triplet, target_size, interpolation=cv2.INTER_LINEAR)\n        processed_24.append(resized_triplet.transpose(2, 0, 1))\n        \n    return np.stack(processed_24, axis=0).astype(np.uint8)\n\n\ndef process_study(study_id):\n    try:\n        # We need the series descriptions for this study from train_series.csv\n        study_df = series_df[series_df['StudyInstanceUID'] == study_id]\n        if len(study_df) == 0: return False\n        \n        out_study_dir = os.path.join(OUT_DIR, str(study_id))\n        os.makedirs(out_study_dir, exist_ok=True)\n        \n        physics_meta = []\n        for _, row in study_df.iterrows():\n            series_id = str(row['SeriesInstanceUID'])\n            plane = str(row['Anatomical_Plane'])\n            fluid_sens = float(row['Fluid_Sensitive'])\n            fat_supp = float(row['Fat_Suppression'])\n            \n            # Map Anatomical_Plane to simple key\n            if plane == 'Sagittal': key = 'Sag'\n            elif plane == 'Coronal': key = 'Cor'\n            elif plane == 'Axial': key = 'Ax'\n            else: continue\n                \n            dicom_dir = os.path.join(DATA_DIR, str(study_id), series_id)\n            if not os.path.exists(dicom_dir): continue\n                \n            dicom_paths = glob.glob(os.path.join(dicom_dir, \"*.dcm\"))\n            \n            # Process and save\n            vol_npy = process_dicom_volume(dicom_paths, laterality=\"R\", target_size=TARGET_SIZE, num_save_slices=NUM_SLICES)\n            \n            # Save as {key}_{series_id}.npy (e.g. Sag_123456.npy)\n            np.save(os.path.join(out_study_dir, f\"{key}_{series_id}.npy\"), vol_npy)\n            \n            # Save physics meta\n            physics_meta.append({\n                'series_id': series_id,\n                'key': key,\n                'fluid_sens': fluid_sens,\n                'fat_supp': fat_supp\n            })\n            \n        pd.DataFrame(physics_meta).to_csv(os.path.join(out_study_dir, \"physics_meta.csv\"), index=False)\n        return True\n    except Exception as e:\n        print(f\"Error on {study_id}: {e}\")\n        return False\n\nif __name__ == \"__main__\":\n    WORKERS = 4\n    \n    # KAGGLE 20GB DISK LIMIT WORKAROUND:\n    # Since the full cache is ~100GB, we must process in chunks of 800 studies per notebook run.\n    # Change these values for each run:\n    # Run 1: 0 to 800\n    # Run 2: 800 to 1600\n    # Run 3: 1600 to 2400\n    # Run 4: 2400 to 3200\n    # Run 5: 3200 to 4000\n    # Run 6: 4000 to 99999\n    START_IDX = 3100\n    END_IDX = 3200\n    \n    os.makedirs(OUT_DIR, exist_ok=True)\n    \n    print(\"Loading train_series.csv...\")\n    series_df = pd.read_csv(DESC_CSV)\n    all_studies = series_df['StudyInstanceUID'].unique().tolist()\n    \n    # Slice the studies to prevent Out Of Space error\n    all_studies = all_studies[START_IDX:END_IDX]\n    \n    print(f\"Found {len(all_studies)} studies for this chunk. Starting multiprocessing with {WORKERS} workers...\")\n    \n    with Pool(WORKERS) as p:\n        results = list(tqdm(p.imap(process_study, all_studies), total=len(all_studies)))\n        \n    print(f\"Successfully processed {sum(results)} / {len(all_studies)} studies.\")\n    print(f\"Cache chunk saved to {OUT_DIR}. Please save this Kaggle output as a Dataset!\")\n\n"}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python"}},"nbformat":4,"nbformat_minor":4}