{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.11","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":91249,"databundleVersionId":11294684,"sourceType":"competition"}],"dockerImageVersionId":31040,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"%load_ext autoreload\n%autoreload 2","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-19T09:00:55.235233Z","iopub.execute_input":"2025-05-19T09:00:55.235756Z","iopub.status.idle":"2025-05-19T09:00:55.267620Z","shell.execute_reply.started":"2025-05-19T09:00:55.235723Z","shell.execute_reply":"2025-05-19T09:00:55.266696Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import polars as pl\n\nlabel_df = pl.read_csv(\"/kaggle/input/byu-locating-bacterial-flagellar-motors-2025/train_labels.csv\")\nlabel_df = label_df.with_columns(\n    volume=pl.col(\"Array shape (axis 0)\") * pl.col(\"Array shape (axis 1)\") * pl.col(\"Array shape (axis 2)\")\n)\ntomo_ids = label_df[\"tomo_id\"].unique(maintain_order=True).to_numpy()","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-05-19T09:00:55.573409Z","iopub.execute_input":"2025-05-19T09:00:55.574379Z","iopub.status.idle":"2025-05-19T09:00:56.780196Z","shell.execute_reply.started":"2025-05-19T09:00:55.574343Z","shell.execute_reply":"2025-05-19T09:00:56.779372Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"label_df","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-19T07:54:59.112480Z","iopub.execute_input":"2025-05-19T07:54:59.112846Z","iopub.status.idle":"2025-05-19T07:54:59.138776Z","shell.execute_reply.started":"2025-05-19T07:54:59.112822Z","shell.execute_reply":"2025-05-19T07:54:59.137085Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from pathlib import Path\n\n\nraw_tomo_path = Path(\"/kaggle/input/byu-locating-bacterial-flagellar-motors-2025/train\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-19T07:51:51.274593Z","iopub.execute_input":"2025-05-19T07:51:51.274952Z","iopub.status.idle":"2025-05-19T07:51:51.295047Z","shell.execute_reply.started":"2025-05-19T07:51:51.274929Z","shell.execute_reply":"2025-05-19T07:51:51.294085Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from pathlib import Path\nimport numpy as np\nimport imageio.v3 as imageio\n\n\ndef load_tomo(tomo_id, tomo_dir, stride=1):\n    tomo_path = Path(tomo_dir) / tomo_id\n    # Load the tomogram\n    voxel = []\n    for path in sorted(tomo_path.glob(\"slice_*.jpg\")):\n        slice_id = int(path.stem.split(\"_\")[-1])\n        if slice_id % stride != 0:\n            continue\n        try:\n            tomogram = np.array(imageio.imread(path))\n            voxel.append(tomogram)\n        except FileNotFoundError:\n            continue\n    voxel = np.stack(voxel, axis=0)\n    return voxel","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-19T07:51:51.461675Z","iopub.execute_input":"2025-05-19T07:51:51.462046Z","iopub.status.idle":"2025-05-19T07:51:51.485055Z","shell.execute_reply.started":"2025-05-19T07:51:51.462024Z","shell.execute_reply":"2025-05-19T07:51:51.482931Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from tqdm import tqdm\n\n\ndef check_outplier(raw_tomo_path, tomo_ids, stride=50):\n    outlier_lefts = []\n    outlier_rights = []\n    outlier_scores = []\n    for tomo_id in tqdm(tomo_ids):\n        voxel = load_tomo(\n            tomo_id,\n            raw_tomo_path,\n            stride=stride,\n        )  # Load the first tomogram\n        count, left = np.histogram(voxel.flatten(), bins=255, range=(0, 255))\n        outlier_lefts.append(count[:50].sum() / count.sum())\n        outlier_rights.append(count[-50:].sum() / count.sum())\n        outlier_scores.append((count[:50].sum() + count[-50].sum()) / count.sum())\n    \n    outlier_lefts = np.array(outlier_lefts)\n    outlier_rights = np.array(outlier_rights)\n    outlier_scores = np.array(outlier_scores)\n    return outlier_scores","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-19T07:51:51.620329Z","iopub.execute_input":"2025-05-19T07:51:51.620695Z","iopub.status.idle":"2025-05-19T07:51:51.640616Z","shell.execute_reply.started":"2025-05-19T07:51:51.620672Z","shell.execute_reply":"2025-05-19T07:51:51.639533Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pickle\n\n\ndef load_pickle(pickle_path):\n    with open(pickle_path, \"rb\") as f:\n        return pickle.load(f)\n\ndef save_pickle(data, pickle_path):\n    with open(pickle_path, \"wb\") as f:\n        pickle.dump(data, f, protocol=pickle.HIGHEST_PROTOCOL)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-19T07:51:54.719123Z","iopub.execute_input":"2025-05-19T07:51:54.719425Z","iopub.status.idle":"2025-05-19T07:51:54.738551Z","shell.execute_reply.started":"2025-05-19T07:51:54.719408Z","shell.execute_reply":"2025-05-19T07:51:54.737529Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"outlier_path = Path(\"outlier_scores.pkl\")\noutlier_scores = check_outplier(raw_tomo_path, tomo_ids, stride=50)\nsave_pickle(outlier_scores, outlier_path)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-19T08:08:40.985693Z","iopub.execute_input":"2025-05-19T08:08:40.986004Z","iopub.status.idle":"2025-05-19T08:08:41.010639Z","shell.execute_reply.started":"2025-05-19T08:08:40.985984Z","shell.execute_reply":"2025-05-19T08:08:41.009620Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import matplotlib.pyplot as plt\n\nplt.hist(outlier_scores, bins=100)\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-19T08:08:48.610269Z","iopub.execute_input":"2025-05-19T08:08:48.610625Z","iopub.status.idle":"2025-05-19T08:08:48.880315Z","shell.execute_reply.started":"2025-05-19T08:08:48.610602Z","shell.execute_reply":"2025-05-19T08:08:48.879349Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"THRESHOLD = 0.5\n\n\nanomaly_idxs = np.where(outlier_scores > THRESHOLD)[0]\nanomaly_tomo_ids = tomo_ids[anomaly_idxs]\nprint(anomaly_tomo_ids)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-19T07:53:04.683615Z","iopub.execute_input":"2025-05-19T07:53:04.684153Z","iopub.status.idle":"2025-05-19T07:53:04.713052Z","shell.execute_reply.started":"2025-05-19T07:53:04.684118Z","shell.execute_reply":"2025-05-19T07:53:04.711221Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"save_pickle(anomaly_tomo_ids, \"anomaly_tomo_ids.pkl\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-19T07:53:04.938407Z","iopub.execute_input":"2025-05-19T07:53:04.938772Z","iopub.status.idle":"2025-05-19T07:53:04.966963Z","shell.execute_reply.started":"2025-05-19T07:53:04.938748Z","shell.execute_reply":"2025-05-19T07:53:04.965769Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"volume_df = label_df.filter(pl.col(\"tomo_id\").is_in(anomaly_tomo_ids)).group_by(\"tomo_id\").agg(pl.all().first())\ntotal_volume = volume_df[\"volume\"].sum()\nprint(f\"{total_volume=}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-19T09:01:19.895440Z","iopub.execute_input":"2025-05-19T09:01:19.895745Z","iopub.status.idle":"2025-05-19T09:01:19.982949Z","shell.execute_reply.started":"2025-05-19T09:01:19.895721Z","shell.execute_reply":"2025-05-19T09:01:19.981918Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"volume_df.describe()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-19T09:01:39.336197Z","iopub.execute_input":"2025-05-19T09:01:39.336529Z","iopub.status.idle":"2025-05-19T09:01:39.355863Z","shell.execute_reply.started":"2025-05-19T09:01:39.336502Z","shell.execute_reply":"2025-05-19T09:01:39.354835Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def percentile_norm(image, lb=1, ub=99):\n    \"\"\"\n    Normalize the image using percentile normalization.\n    \"\"\"\n    pl, pu = np.percentile(image, (lb, ub))\n    image = (image - pl) / (pu - pl)\n    image = np.clip(image, 0, 1)\n    return image","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-19T07:56:11.098123Z","iopub.execute_input":"2025-05-19T07:56:11.098480Z","iopub.status.idle":"2025-05-19T07:56:11.122212Z","shell.execute_reply.started":"2025-05-19T07:56:11.098457Z","shell.execute_reply":"2025-05-19T07:56:11.120548Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from pathlib import Path\nimport numpy as np\nimport imageio.v3 as imageio\n\n\ndef load_slice(tomo_id, tomo_dir, slice_id):\n    slice_path = Path(tomo_dir) / tomo_id / f\"slice_{slice_id:04d}.jpg\"\n    image = np.array(imageio.imread(slice_path))\n    return image","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-19T07:56:11.278889Z","iopub.execute_input":"2025-05-19T07:56:11.279209Z","iopub.status.idle":"2025-05-19T07:56:11.300139Z","shell.execute_reply.started":"2025-05-19T07:56:11.279188Z","shell.execute_reply":"2025-05-19T07:56:11.298753Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def correct_tomo(x):\n    x = (x.astype(np.uint8) + 127).astype(np.float32)\n    return x","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-19T07:56:11.469625Z","iopub.execute_input":"2025-05-19T07:56:11.470689Z","iopub.status.idle":"2025-05-19T07:56:11.491533Z","shell.execute_reply.started":"2025-05-19T07:56:11.470634Z","shell.execute_reply":"2025-05-19T07:56:11.490509Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"for tomo_id in anomaly_tomo_ids:\n    tomo_path = raw_tomo_path / tomo_id\n    depth = len(list(tomo_path.glob(\"*.jpg\")))\n    \n    slice_id = depth // 2\n    image = load_slice(tomo_id, raw_tomo_path, slice_id)\n    image_org = image.copy()\n    image = correct_tomo(image)\n\n    image_norm = (percentile_norm(image, 1, 99) * 255).astype(np.uint8)\n\n    print(f\"{tomo_id=}\")\n\n    _, ax = plt.subplots()\n    ax.hist(image_org.flatten(), bins=255, range=(0, 255))\n    plt.show()\n\n    _, ax = plt.subplots()\n    ax.hist(image.flatten(), bins=255, range=(0, 255))\n    plt.show()\n\n    _, ax = plt.subplots()\n    ax.hist(image_norm.flatten(), bins=255)\n    plt.show()\n\n    _, ax = plt.subplots()\n    ax.imshow(image_org, cmap=\"gray\")\n    plt.show()\n\n    _, ax = plt.subplots()\n    ax.imshow(image, cmap=\"gray\")\n    plt.show()\n\n    _, ax = plt.subplots()\n    ax.imshow(image_norm, cmap=\"gray\")\n    plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-19T07:56:11.623502Z","iopub.execute_input":"2025-05-19T07:56:11.624623Z","iopub.status.idle":"2025-05-19T07:56:50.289774Z","shell.execute_reply.started":"2025-05-19T07:56:11.624582Z","shell.execute_reply":"2025-05-19T07:56:50.288795Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def calc_hist(x_hist, lb=1, ub=99, bins=512):\n    hist, bins = np.histogram(x_hist, bins=bins)\n    cdf = np.cumsum(hist) / hist.sum()\n    l_idx = np.searchsorted(cdf, lb / 100)\n    u_idx = np.searchsorted(cdf, ub / 100)\n    lower, upper = bins[l_idx], bins[u_idx]\n    return lower, upper\n\n\ndef process_volume(x, lower, upper):\n    x = correct_tomo(x)\n    \n    x = np.clip(x, lower, upper)\n    x = (x - lower) / (upper - lower)\n\n    # Quantize\n    x = np.round(x.clip(0, 1) * 255).astype(np.uint8)\n\n    return x","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-19T07:44:16.292025Z","iopub.execute_input":"2025-05-19T07:44:16.292409Z","iopub.status.idle":"2025-05-19T07:44:16.314478Z","shell.execute_reply.started":"2025-05-19T07:44:16.292387Z","shell.execute_reply":"2025-05-19T07:44:16.313269Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"out_dir = Path(\"volumes\")\nout_dir.mkdir(exist_ok=True, parents=True)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-19T07:44:17.866012Z","iopub.execute_input":"2025-05-19T07:44:17.866306Z","iopub.status.idle":"2025-05-19T07:44:17.884922Z","shell.execute_reply.started":"2025-05-19T07:44:17.866285Z","shell.execute_reply":"2025-05-19T07:44:17.883445Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\n\nDEBUG = os.environ.get(\"KAGGLE_KERNEL_RUN_TYPE\") == 'Interactive'\nDEBUG","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-19T07:44:18.109137Z","iopub.execute_input":"2025-05-19T07:44:18.109422Z","iopub.status.idle":"2025-05-19T07:44:18.129374Z","shell.execute_reply.started":"2025-05-19T07:44:18.109401Z","shell.execute_reply":"2025-05-19T07:44:18.128091Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from tqdm.auto import tqdm\nfrom PIL import Image\n\n\ntarget_tomo_ids = anomaly_tomo_ids\nif DEBUG:\n    target_tomo_ids = [\"tomo_3b8291\", \"tomo_b18127\"]\n    target_tomo_ids = [\"tomo_b18127\"]\nfor tomo_id in tqdm(target_tomo_ids):\n    voxel = load_tomo(\n        tomo_id,\n        raw_tomo_path,\n        stride=1,\n    )\n    lower, upper = calc_hist(voxel.flatten(), lb=1, ub=99, bins=512)\n    voxel_norm = process_volume(voxel, lower, upper)\n\n    voxel_path = out_dir / tomo_id\n    voxel_path.mkdir(exist_ok=True, parents=True)\n    for slice_id in tqdm(range(voxel.shape[0])):\n        slice_path = voxel_path / f\"slice_{slice_id:04d}.jpg\"\n        slice_path.parent.mkdir(exist_ok=True, parents=True)\n        img = Image.fromarray(voxel_norm[slice_id])\n        img.save(slice_path, format=\"JPEG\", quality=50)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-19T07:46:40.729000Z","iopub.execute_input":"2025-05-19T07:46:40.729883Z","iopub.status.idle":"2025-05-19T07:47:03.234453Z","shell.execute_reply.started":"2025-05-19T07:46:40.729818Z","shell.execute_reply":"2025-05-19T07:47:03.233404Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from matplotlib.patches import Rectangle\n\n\nfor tomo_id in anomaly_tomo_ids:\n    df = label_df.filter(pl.col(\"tomo_id\") == tomo_id)\n    for row in df.iter_rows(named=True):\n        has_motor = row[\"Number of motors\"] > 0\n        z = (\n            row[\"Motor axis 0\"] if has_motor else row[\"Array shape (axis 0)\"] // 2\n        )\n        y = row[\"Motor axis 1\"] if has_motor else 0\n        x = row[\"Motor axis 2\"] if has_motor else 0\n        voxel_spacing = row[\"Voxel spacing\"]\n\n        try:\n            image = load_slice(\n                tomo_id,\n                out_dir,\n                int(z),\n            )\n            ax = plt.gca()\n            ax.hist(image.flatten(), bins=255)\n            ax.set(title=f\"{tomo_id=}, {z=}, {y=}, {x=}\")\n            plt.show()\n    \n            fig, ax = plt.subplots(figsize=(8, 8))\n            ax.imshow(image, cmap=\"gray\")\n            if has_motor:\n                s = 1000 / voxel_spacing\n                ax.add_patch(\n                    Rectangle(\n                        (x - s // 2, y - s // 2),\n                        s,\n                        s,\n                        linewidth=0.5,\n                        edgecolor=\"r\",\n                        facecolor=\"none\",\n                    )\n                )\n                ax.set(title=f\"{tomo_id=}, {z=}, {y=}, {x=}\")\n            plt.show()\n        except Exception as e:\n            print(f\"⚠️ Failed to process {tomo_id}: {e}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-19T08:00:03.550473Z","iopub.execute_input":"2025-05-19T08:00:03.550889Z","iopub.status.idle":"2025-05-19T08:00:04.631083Z","shell.execute_reply.started":"2025-05-19T08:00:03.550859Z","shell.execute_reply":"2025-05-19T08:00:04.629765Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!du -sh","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-19T08:00:04.633100Z","iopub.execute_input":"2025-05-19T08:00:04.633500Z","iopub.status.idle":"2025-05-19T08:00:04.784597Z","shell.execute_reply.started":"2025-05-19T08:00:04.633475Z","shell.execute_reply":"2025-05-19T08:00:04.783115Z"}},"outputs":[],"execution_count":null}]}