{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Load 🤖 convert DICOM to 3D volume","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"!pip download -q \"python-gdcm\" pydicom pylibjpeg \"opencv-python-headless\" --dest frozen_packages --prefer-binary\n!pip wheel -q https://github.com/Borda/kaggle_vol-3D-classify/archive/refs/heads/main.zip --wheel-dir frozen_packages --prefer-binary\n!rm frozen_packages/torch-*\n!ls -lh frozen_packages","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-01T21:49:46.859345Z","iopub.execute_input":"2022-08-01T21:49:46.859806Z","iopub.status.idle":"2022-08-01T21:50:33.858683Z","shell.execute_reply.started":"2022-08-01T21:49:46.859771Z","shell.execute_reply":"2022-08-01T21:50:33.857063Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install -qU \"python-gdcm\" pydicom pylibjpeg kaggle_vol3d_classify --find-links frozen_packages --no-index","metadata":{"execution":{"iopub.status.busy":"2022-08-01T21:50:33.861779Z","iopub.execute_input":"2022-08-01T21:50:33.862235Z","iopub.status.idle":"2022-08-01T21:50:46.296573Z","shell.execute_reply.started":"2022-08-01T21:50:33.862178Z","shell.execute_reply":"2022-08-01T21:50:46.294847Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%matplotlib inline\n\nimport os\nimport glob\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nPATH_DATASET = \"/kaggle/input/rsna-2022-cervical-spine-fracture-detection\"","metadata":{"execution":{"iopub.status.busy":"2022-08-01T21:23:45.509633Z","iopub.execute_input":"2022-08-01T21:23:45.51054Z","iopub.status.idle":"2022-08-01T21:23:45.521265Z","shell.execute_reply.started":"2022-08-01T21:23:45.510495Z","shell.execute_reply":"2022-08-01T21:23:45.519552Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Loading & saving DICOM image","metadata":{}},{"cell_type":"code","source":"import cv2\nimport pydicom\nimport torch\nfrom PIL import Image\nfrom dipy.denoise.nlmeans import nlmeans\nfrom dipy.denoise.noise_estimate import estimate_sigma\nfrom pydicom.pixel_data_handlers import apply_voi_lut\nfrom kaggle_volclassif.utils import interpolate_volume\nfrom skimage import exposure\n    \n\ndef convert_volume(dir_path: str, out_dir: str = \"train_volumes\", size = (224, 224, 224)):\n    ls_imgs = glob.glob(os.path.join(dir_path, \"*.dcm\"))\n    ls_imgs = sorted(ls_imgs, key=lambda p: int(os.path.splitext(os.path.basename(p))[0]))\n\n    imgs = []\n    for p_img in ls_imgs:\n        dicom = pydicom.dcmread(p_img)\n        img = apply_voi_lut(dicom.pixel_array, dicom)\n        img = cv2.resize(img, size[:2], interpolation=cv2.INTER_LINEAR)\n        imgs.append(img.tolist())\n    vol = torch.tensor(imgs, dtype=torch.float32)\n\n    vol = (vol - vol.min()) / float(vol.max() - vol.min())\n    vol = interpolate_volume(vol, size).numpy()\n    \n    # https://scikit-image.org/docs/stable/auto_examples/color_exposure/plot_adapt_hist_eq_3d.html\n    vol = exposure.equalize_adapthist(vol, kernel_size=np.array([64, 64, 64]), clip_limit=0.01)\n    # vol = exposure.equalize_hist(vol)\n    vol = np.clip(vol * 255, 0, 255).astype(np.uint8)\n    \n    path_npz = os.path.join(out_dir, f\"{os.path.basename(dir_path)}.npz\")\n    np.savez_compressed(path_npz, vol)","metadata":{"execution":{"iopub.status.busy":"2022-08-01T22:10:27.34954Z","iopub.execute_input":"2022-08-01T22:10:27.350922Z","iopub.status.idle":"2022-08-01T22:10:27.366033Z","shell.execute_reply.started":"2022-08-01T22:10:27.350871Z","shell.execute_reply":"2022-08-01T22:10:27.364767Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Process all images 🤖","metadata":{}},{"cell_type":"code","source":"from pprint import pprint\nfrom joblib import Parallel, delayed\nfrom tqdm.auto import tqdm\n\n! rm -rf train_volumes\n! mkdir train_volumes\n\nls_dirs = [p for p in glob.glob(os.path.join(PATH_DATASET, \"train_images\", \"*\")) if os.path.isdir(p)]\nprint(f\"volumes: {len(ls_dirs)}\")\n\n_= Parallel(n_jobs=4)(delayed(convert_volume)(p_dir) for p_dir in tqdm(ls_dirs))\n\n! ls -lh train_volumes","metadata":{"execution":{"iopub.status.busy":"2022-08-01T22:10:29.31101Z","iopub.execute_input":"2022-08-01T22:10:29.312134Z"},"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Show 🔎 few samples","metadata":{}},{"cell_type":"code","source":"from ipywidgets import interact, IntSlider\nfrom matplotlib.patches import PathPatch, Rectangle\nfrom matplotlib.path import Path\n\n\ndef _draw_line(ax, coords, clr='g'):\n    line = Path(coords, [Path.MOVETO, Path.LINETO])\n    pp = PathPatch(line, linewidth=3, edgecolor=clr, facecolor='none')\n    ax.add_patch(pp)\n\ndef _set_axes_labels(ax, axes_x, axes_y):\n    ax.set_xlabel(axes_x)\n    ax.set_ylabel(axes_y)\n    ax.set_aspect('equal', 'box')\n\n_rec_prop = dict(linewidth=5, facecolor='none')\n\ndef show_volume(vol, z, y, x, fig_size=(9, 9)):\n    fig, axarr = plt.subplots(nrows=2, ncols=2, figsize=fig_size)\n    v_z, v_y, v_x = vol.shape\n    axarr[0, 0].imshow(vol[z, :, :], cmap=\"gray\")\n    axarr[0, 0].add_patch(Rectangle((-1, -1), v_x, v_y, edgecolor='r', **_rec_prop))\n    _draw_line(axarr[0, 0], [(x, 0), (x, v_y)], \"g\")\n    _draw_line(axarr[0, 0], [(0, y), (v_x, y)], \"b\")\n    _set_axes_labels(axarr[0, 0], \"X\", \"Y\")\n    axarr[0, 1].imshow(vol[:, :, x].T, cmap=\"gray\")\n    axarr[0, 1].add_patch(Rectangle((-1, -1), v_z, v_y, edgecolor='g', **_rec_prop))\n    _draw_line(axarr[0, 1], [(z, 0), (z, v_y)], \"r\")\n    _draw_line(axarr[0, 1], [(0, y), (v_x, y)], \"b\")\n    _set_axes_labels(axarr[0, 1], \"Z\", \"Y\")\n    im = axarr[1, 0].imshow(vol[:, y, :], cmap=\"gray\")\n    axarr[1, 0].add_patch(Rectangle((-1, -1), v_x, v_z, edgecolor='b', **_rec_prop))\n    _draw_line(axarr[1, 0], [(0, z), (v_x, z)], \"r\")\n    _draw_line(axarr[1, 0], [(x, 0), (x, v_y)], \"g\")\n    _set_axes_labels(axarr[1, 0], \"X\", \"Z\")\n    plt.colorbar(im, ax=axarr[1, 1])\n    axarr[1, 1].set_axis_off()\n    fig.tight_layout()\n\n\ndef interactive_show(volume):\n    vol_shape = volume.shape\n    interact(\n        lambda x, y, z: plt.show(show_volume(volume, z, y, x)),\n        z=IntSlider(min=0, max=vol_shape[0], step=2, value=int(vol_shape[0] / 2)),\n        y=IntSlider(min=0, max=vol_shape[1], step=5, value=int(vol_shape[1] / 2)),\n        x=IntSlider(min=0, max=vol_shape[2], step=5, value=int(vol_shape[2] / 2)),\n    )","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-01T22:08:55.392573Z","iopub.execute_input":"2022-08-01T22:08:55.393105Z","iopub.status.idle":"2022-08-01T22:08:55.420717Z","shell.execute_reply.started":"2022-08-01T22:08:55.393065Z","shell.execute_reply":"2022-08-01T22:08:55.419475Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ls_vols = glob.glob(os.path.join(\"train_volumes\", \"*.npz\"))\n\nvol = np.load(ls_vols[0])['arr_0']\n\ninteractive_show(vol)","metadata":{"execution":{"iopub.status.busy":"2022-08-01T22:08:58.693264Z","iopub.execute_input":"2022-08-01T22:08:58.694468Z","iopub.status.idle":"2022-08-01T22:08:59.641346Z","shell.execute_reply.started":"2022-08-01T22:08:58.69438Z","shell.execute_reply":"2022-08-01T22:08:59.640199Z"},"trusted":true},"execution_count":null,"outputs":[]}]}