{"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-09-06T01:23:48.728192Z","iopub.execute_input":"2022-09-06T01:23:48.73007Z","iopub.status.idle":"2022-09-06T01:24:21.446948Z","shell.execute_reply.started":"2022-09-06T01:23:48.729976Z","shell.execute_reply":"2022-09-06T01:24:21.443624Z"},"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-09-06T01:24:21.448879Z","iopub.status.idle":"2022-09-06T01:24:21.450337Z","shell.execute_reply.started":"2022-09-06T01:24:21.45008Z","shell.execute_reply":"2022-09-06T01:24:21.450111Z"},"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-09-06T01:24:21.45166Z","iopub.status.idle":"2022-09-06T01:24:21.452369Z","shell.execute_reply.started":"2022-09-06T01:24:21.452138Z","shell.execute_reply":"2022-09-06T01:24:21.452162Z"},"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 = \"test_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    z_and_imgs = []\n    for p_img in ls_imgs:\n        dicom = pydicom.dcmread(p_img)\n        b = dicom.RescaleIntercept\n        m = dicom.RescaleSlope\n        img = m * dicom.pixel_array + b\n        img = cv2.resize(img, size[:2], interpolation=cv2.INTER_LINEAR)\n        z_pos = float(dicom.ImagePositionPatient[-1])\n        z_and_imgs.append([z_pos,img])\n    \n    z_and_imgs.sort(key=lambda tup: tup[0])  # sort using z as index\n    n = 1 # tuple index of img\n    imgs = [x[n].tolist() for x in z_and_imgs]\n    vol = torch.tensor(imgs, dtype=torch.float32)\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 = vol.astype(np.float16)\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-09-06T01:24:21.453755Z","iopub.status.idle":"2022-09-06T01:24:21.454244Z","shell.execute_reply.started":"2022-09-06T01:24:21.454016Z","shell.execute_reply":"2022-09-06T01:24:21.454054Z"},"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! rm -rf test_volumes\n! mkdir test_volumes\n\nls_dirs = [p for p in glob.glob(os.path.join(PATH_DATASET, \"test_images\", \"*\")) if os.path.isdir(p)]\nls_dirs.sort()\n#ls_dirs = ls_dirs[:]\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 test_volumes","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-09-06T01:24:21.456345Z","iopub.status.idle":"2022-09-06T01:24:21.456851Z","shell.execute_reply.started":"2022-09-06T01:24:21.456638Z","shell.execute_reply":"2022-09-06T01:24:21.45666Z"},"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-09-06T01:24:21.459227Z","iopub.status.idle":"2022-09-06T01:24:21.459715Z","shell.execute_reply.started":"2022-09-06T01:24:21.459492Z","shell.execute_reply":"2022-09-06T01:24:21.459519Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ls_vols = glob.glob(os.path.join(\"test_volumes\", \"*.npz\"))\n\nvol = np.load(ls_vols[0])['arr_0']\n\ninteractive_show(vol)","metadata":{"execution":{"iopub.status.busy":"2022-09-06T01:24:21.461657Z","iopub.status.idle":"2022-09-06T01:24:21.462198Z","shell.execute_reply.started":"2022-09-06T01:24:21.461957Z","shell.execute_reply":"2022-09-06T01:24:21.461985Z"},"trusted":true},"execution_count":null,"outputs":[]}]}