{"metadata":{"accelerator":"GPU","colab":{"authorship_tag":"ABX9TyMVWrwOAJwKc4ndWcRnVD0b","gpuType":"V100","machine_shape":"hm","provenance":[]},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceId":52254,"databundleVersionId":6863140,"sourceType":"competition"},{"sourceId":154489631,"sourceType":"kernelVersion"}],"dockerImageVersionId":30615,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"!pip install -q opencv-python\n!pip install -q dicomsdl\n!pip install -q -U diffdrr\n!pip install torch==1.13.1+cu116 torchvision==0.14.1+cu116 torchaudio==0.13.1 --extra-index-url https://download.pytorch.org/whl/cu116","metadata":{"executionInfo":{"elapsed":61514,"status":"ok","timestamp":1699516365546,"user":{"displayName":"competition competition","userId":"10960398843493506800"},"user_tz":-540},"id":"J0PU-Zm57jWi","outputId":"8c0544e0-c830-420e-b051-97f76bf1c869","execution":{"iopub.status.busy":"2024-01-15T07:23:11.404291Z","iopub.execute_input":"2024-01-15T07:23:11.404588Z","iopub.status.idle":"2024-01-15T07:25:51.872848Z","shell.execute_reply.started":"2024-01-15T07:23:11.404563Z","shell.execute_reply":"2024-01-15T07:25:51.871558Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import sys\nimport torch\npyt_version_str=torch.__version__.split(\"+\")[0].replace(\".\", \"\")\nversion_str=\"\".join([\n    f\"py3{sys.version_info.minor}_cu\",\n    torch.version.cuda.replace(\".\",\"\"),\n    f\"_pyt{pyt_version_str}\"\n])\n!pip install fvcore iopath\n!pip install --no-index --no-cache-dir pytorch3d -f https://dl.fbaipublicfiles.com/pytorch3d/packaging/wheels/{version_str}/download.html","metadata":{"executionInfo":{"elapsed":16520,"status":"ok","timestamp":1699516382061,"user":{"displayName":"competition competition","userId":"10960398843493506800"},"user_tz":-540},"id":"SVOjkGfd7jYp","outputId":"a3c43500-1b69-4772-f827-b3d1c6968b8b","execution":{"iopub.status.busy":"2024-01-15T07:25:51.874975Z","iopub.execute_input":"2024-01-15T07:25:51.875303Z","iopub.status.idle":"2024-01-15T07:26:27.58484Z","shell.execute_reply.started":"2024-01-15T07:25:51.875273Z","shell.execute_reply":"2024-01-15T07:26:27.583752Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import tensorflow as tf\nimport zipfile\nimport os\nimport torch\nimport pydicom\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport cv2\nimport glob\nimport dicomsdl\nimport matplotlib.pyplot as plt\nimport torch\nfrom tqdm import tqdm\nimport gc\nfrom PIL import Image\nimport nibabel as nib\n\nfrom diffdrr.drr import DRR\nfrom diffdrr.data import load_example_ct\n# from diffdrr.visualization import plot_drr","metadata":{"executionInfo":{"elapsed":4249,"status":"ok","timestamp":1699516386303,"user":{"displayName":"competition competition","userId":"10960398843493506800"},"user_tz":-540},"id":"YDsateFJ7jca","outputId":"fddfa62c-047c-4907-e836-7300028e1610","execution":{"iopub.status.busy":"2024-01-15T07:26:27.586913Z","iopub.execute_input":"2024-01-15T07:26:27.587681Z","iopub.status.idle":"2024-01-15T07:26:41.341764Z","shell.execute_reply.started":"2024-01-15T07:26:27.587641Z","shell.execute_reply":"2024-01-15T07:26:41.340974Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import shutil\nfrom scipy import ndimage\nfrom scipy.ndimage import zoom\nfrom ipywidgets import IntSlider, interact\nfrom matplotlib import animation, rc\nfrom matplotlib.patches import PathPatch, Rectangle\nfrom matplotlib.path import Path\nimport cv2\nfrom IPython.display import HTML\nimport math\nimport pandas as pd\nimport random","metadata":{"execution":{"iopub.status.busy":"2024-01-15T07:26:41.343968Z","iopub.execute_input":"2024-01-15T07:26:41.344515Z","iopub.status.idle":"2024-01-15T07:26:41.354186Z","shell.execute_reply.started":"2024-01-15T07:26:41.344487Z","shell.execute_reply":"2024-01-15T07:26:41.353026Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def load_scan_dcm(directory):\n    # DICOMファイルのリストを取得\n    dicom_files = [pydicom.dcmread(os.path.join(directory, f)) for f in os.listdir(directory) if f.endswith('.dcm')]\n    # スライス位置に基づいてファイルをソート\n    dicom_files.sort(key=lambda x: float(x.ImagePositionPatient[2]))\n\n    # スライス間の距離を計算してスライス厚を決定する\n    try:\n        slice_thickness = np.abs(dicom_files[0].ImagePositionPatient[2] - dicom_files[1].ImagePositionPatient[2])\n    except AttributeError:\n        slice_thickness = np.abs(dicom_files[0].SliceLocation - dicom_files[1].SliceLocation)\n\n    # 画像データを格納するための配列を作成\n    image_array = np.stack([file.pixel_array for file in dicom_files])\n    image_array = image_array.astype(np.float32)  # dtypeを変更する場合\n    # img_show(image_array)\n    image_array = np.swapaxes(image_array, 1, 2)\n    image_array = np.flip(image_array, axis=1)\n    # 上下（y軸方向）を反転\n    image_array = np.flip(image_array, axis=0)\n    # img_show(image_array)\n    image_array = np.transpose(image_array, (2, 1, 0))\n\n    # ウィンドウレベルの調整が必要な場合\n    # image_array += 2000\n\n    # スペーシング情報を取得a\n    pixel_spacing = list(dicom_files[0].PixelSpacing)\n    pixel_spacing.append(slice_thickness)\n    spacing = np.array(pixel_spacing, dtype=np.float32)\n\n    return image_array, spacing","metadata":{"executionInfo":{"elapsed":13,"status":"ok","timestamp":1699516386304,"user":{"displayName":"competition competition","userId":"10960398843493506800"},"user_tz":-540},"id":"P9mcqxCx7jdh","execution":{"iopub.status.busy":"2024-01-15T07:26:41.355808Z","iopub.execute_input":"2024-01-15T07:26:41.356237Z","iopub.status.idle":"2024-01-15T07:26:41.380679Z","shell.execute_reply.started":"2024-01-15T07:26:41.356195Z","shell.execute_reply":"2024-01-15T07:26:41.379902Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def load_nii_gz_file(file_path):\n    # NIfTIファイルの読み込み\n    nii_image = nib.load(file_path)\n\n    # 画像データを取得\n    image_array = nii_image.get_fdata()\n    image_array = image_array.astype(np.float32)  # dtypeを変更する場合\n\n    # NIfTIファイルからスペーシング情報を取得\n    header = nii_image.header\n    zooms = header.get_zooms()\n    spacing = np.array(zooms, dtype=np.float32)\n\n    # 必要に応じて画像の向きを調整\n    image_array = np.transpose(image_array, (1, 2, 0))\n    image_array = np.swapaxes(image_array, 1, 2)\n    image_array = np.flip(image_array, axis=2)\n    image_array = np.flip(image_array, axis=0)\n    # image_array = np.transpose(image_array, (2, 1, 0))\n\n    return image_array, spacing","metadata":{"execution":{"iopub.status.busy":"2024-01-15T07:26:41.381904Z","iopub.execute_input":"2024-01-15T07:26:41.382468Z","iopub.status.idle":"2024-01-15T07:26:41.397433Z","shell.execute_reply.started":"2024-01-15T07:26:41.382417Z","shell.execute_reply":"2024-01-15T07:26:41.396618Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_bounding_boxes(heatmap, threshold=0.15, otsu=False):\n    \"\"\"Get bounding boxes from heatmap\"\"\"\n    p_heatmap = np.copy(heatmap)\n    if otsu:\n        # Otsu's thresholding method to find the bounding boxes\n        threshold, p_heatmap = cv2.threshold(\n            heatmap, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU\n        )\n    else:\n        # Using a fixed threshold\n        p_heatmap[p_heatmap < threshold * 255] = 0\n        p_heatmap[p_heatmap >= threshold * 255] = 1\n    # find the contours in the thresholded heatmap\n    contours = cv2.findContours(p_heatmap, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)\n    contours = contours[0] if len(contours) == 2 else contours[1]\n    # get the bounding boxes from the contours\n    bboxes = []\n    for c in contours:\n        x, y, w, h = cv2.boundingRect(c)\n        bboxes.append([x, y, x + w, y + h])\n    return bboxes\ndef get_bbox_patches(bboxes, color='r', linewidth=2):\n    \"\"\"Get patches for bounding boxes\"\"\"\n    patches = []\n    for bbox in bboxes:\n        x1, y1, x2, y2 = bbox\n        patches.append(\n            Rectangle(\n                (x1, y1),\n                x2 - x1,\n                y2 - y1,\n                edgecolor=color,\n                facecolor='none',\n                linewidth=linewidth,\n            )\n        )\n    return patches","metadata":{"execution":{"iopub.status.busy":"2024-01-15T07:26:41.398399Z","iopub.execute_input":"2024-01-15T07:26:41.398697Z","iopub.status.idle":"2024-01-15T07:26:41.411608Z","shell.execute_reply.started":"2024-01-15T07:26:41.398672Z","shell.execute_reply":"2024-01-15T07:26:41.410784Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def create_animation(array, case, heatmap=None, alpha=0.3):\n    \"\"\"Create an animation of a volume\"\"\"\n    # array = np.transpose(array, (2, 0, 1))\n    if heatmap is not None:\n        heatmap = np.transpose(heatmap, (2, 0, 1))\n    fig = plt.figure(figsize=(4, 4))\n    images = []\n    for idx, image in enumerate(array):\n        # plot image without notifying animation\n        image_plot = plt.imshow(image, animated=True, cmap='bone')\n        aux = [image_plot]\n        if heatmap is not None:\n            image_plot2 = plt.imshow(\n                heatmap[idx], animated=True, cmap='jet', alpha=alpha, extent=image_plot.get_extent())\n            aux.append(image_plot2)\n            # add bounding boxes to the heatmap image as animated patches\n            bboxes = get_bounding_boxes(heatmap[idx])\n            patches = get_bbox_patches(bboxes)\n            aux.extend(image_plot2.axes.add_patch(patch) for patch in patches)\n        images.append(aux)\n    plt.axis('off')\n    plt.tight_layout()\n    plt.subplots_adjust(top=0.90)\n    plt.title(f'Patient ID: {case}', fontsize=16)\n    ani = animation.ArtistAnimation(\n        fig, images, interval=5000//len(array), blit=False, repeat_delay=1000)\n    plt.close()\n    return ani","metadata":{"execution":{"iopub.status.busy":"2024-01-15T07:26:41.412718Z","iopub.execute_input":"2024-01-15T07:26:41.413448Z","iopub.status.idle":"2024-01-15T07:26:41.424938Z","shell.execute_reply.started":"2024-01-15T07:26:41.413422Z","shell.execute_reply":"2024-01-15T07:26:41.424239Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(torch.cuda.is_available())","metadata":{"execution":{"iopub.status.busy":"2024-01-15T07:26:41.426092Z","iopub.execute_input":"2024-01-15T07:26:41.426375Z","iopub.status.idle":"2024-01-15T07:26:41.490027Z","shell.execute_reply.started":"2024-01-15T07:26:41.42635Z","shell.execute_reply":"2024-01-15T07:26:41.488982Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mode = 'x-ray'","metadata":{"execution":{"iopub.status.busy":"2024-01-15T07:26:41.492772Z","iopub.execute_input":"2024-01-15T07:26:41.493101Z","iopub.status.idle":"2024-01-15T07:26:41.499628Z","shell.execute_reply.started":"2024-01-15T07:26:41.493074Z","shell.execute_reply":"2024-01-15T07:26:41.498645Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_first_series_paths(base_dir):\n    \"\"\"\n    各patient_idに対して、最初のseries_idのパスを取得します。\n    \n    :param base_dir: 検索を開始するベースディレクトリ\n    :return: 各patient_idの最初のseries_idのパスのリスト\n    \"\"\"\n    first_series_paths = []\n\n    # train_imagesディレクトリ内の各patient_idディレクトリを走査\n    for patient_id in os.listdir(base_dir):\n        patient_path = os.path.join(base_dir, patient_id)\n\n        # patient_idディレクトリがディレクトリであることを確認\n        if os.path.isdir(patient_path):\n            # patient_idディレクトリ内の最初のseries_idディレクトリを見つける\n            for series_id in sorted(os.listdir(patient_path)):\n                series_path = os.path.join(patient_path, series_id)\n                \n                # series_idディレクトリがディレクトリであることを確認\n                if os.path.isdir(series_path):\n                    first_series_paths.append(series_path)\n                    break  # 最初のseries_idのみ必要\n\n    return first_series_paths","metadata":{"execution":{"iopub.status.busy":"2024-01-15T07:26:41.500833Z","iopub.execute_input":"2024-01-15T07:26:41.501105Z","iopub.status.idle":"2024-01-15T07:26:41.510946Z","shell.execute_reply.started":"2024-01-15T07:26:41.501081Z","shell.execute_reply":"2024-01-15T07:26:41.510053Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ベースディレクトリを指定\nbase_directory = '/kaggle/input/rsna-2023-abdominal-trauma-detection/train_images'\n\n# パスを取得\npaths = get_first_series_paths(base_directory)","metadata":{"execution":{"iopub.status.busy":"2024-01-15T07:26:41.512291Z","iopub.execute_input":"2024-01-15T07:26:41.513077Z","iopub.status.idle":"2024-01-15T07:26:55.308946Z","shell.execute_reply.started":"2024-01-15T07:26:41.513048Z","shell.execute_reply":"2024-01-15T07:26:55.307894Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"img_save_dir_xray = \"/kaggle/working/pseudo-x-ray\"\nimg_save_dir_mask = \"/kaggle/working/pseudo-mask-gray\"\n\nif not os.path.exists(img_save_dir_xray):\n    os.makedirs(img_save_dir_xray)\n    \nif not os.path.exists(img_save_dir_mask):\n    os.makedirs(img_save_dir_mask)","metadata":{"execution":{"iopub.status.busy":"2024-01-15T07:26:55.310087Z","iopub.execute_input":"2024-01-15T07:26:55.310364Z","iopub.status.idle":"2024-01-15T07:26:55.317395Z","shell.execute_reply.started":"2024-01-15T07:26:55.31034Z","shell.execute_reply":"2024-01-15T07:26:55.316468Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for path in tqdm(paths[:500]):\n    print(path)\n    input_path_xray = path\n    volume_xray, spacing_xray = load_scan_dcm(input_path_xray)\n    \n    if(volume_xray.shape[0] != 512 or volume_xray.shape[1] != 512 or volume_xray.shape[2] < 400): continue\n\n    input_path_mask = f\"/kaggle/input/totalsegmentator-colon-1/mask/{input_path_xray.split('/')[-2]}\"\n    volume_mask, spacing_mask = load_nii_gz_file(input_path_mask+\"/colon.nii.gz\")\n    \n    bx_xray, by_xray = torch.tensor(volume_xray.shape[:2]) * torch.tensor(spacing_xray[:2]) / 2\n    bz_xray = torch.tensor(volume_xray.shape[2]) * torch.tensor(spacing_xray[2]) * 0.5\n    volume_xray = np.ascontiguousarray(volume_xray)\n\n    bx_mask, by_mask = torch.tensor(volume_mask.shape[:2]) * torch.tensor(spacing_mask[:2]) / 2\n    bz_mask = torch.tensor(volume_mask.shape[2]) * torch.tensor(spacing_mask[2]) * 0.5\n    volume_mask = np.ascontiguousarray(volume_mask)\n\n    # Initialize the DRR module for generating synthetic X-rays\n    device = \"cuda\" if torch.cuda.is_available() else \"cpu\"\n    drr_xray = DRR(\n        volume_xray,      # The CT volume as a numpy array\n        spacing_xray,     # Voxel dimensions of the CT\n        sdr=500,   # Source-to-detector radius (half of the source-to-detector distance)\n        height=int(volume_xray.shape[2]/2),  # Height of the DRR (if width is not seperately provided, the generated image is square)\n        width=int(volume_xray.shape[0]/2),\n        delx=3.0,    # Pixel spacing (in mm)\n    ).to(device)\n\n    drr_mask = DRR(\n        volume_mask,      # The CT volume as a numpy array\n        spacing_mask,     # Voxel dimensions of the CT\n        sdr=500,   # Source-to-detector radius (half of the source-to-detector distance)\n        height=int(volume_mask.shape[2]/2),  # Height of the DRR (if width is not seperately provided, the generated image is square)\n        width=int(volume_mask.shape[0]/2),\n        delx=3.0,    # Pixel spacing (in mm)\n    ).to(device)\n\n     # Set the camera pose with rotation (yaw, pitch, roll) and translation (x, y, z)\n    rotation = torch.tensor([[torch.pi, 0.0, torch.pi / 2]], device=device)\n    translation_xray = torch.tensor([[bx_xray, by_xray, bz_xray]], device=device)\n    translation_mask = torch.tensor([[bx_mask, by_mask, bz_mask]], device=device)\n\n    # 📸 Also note that DiffDRR can take many representations of SO(3) 📸\n    # For example, quaternions, rotation matrix, axis-angle, etc...\n    try:\n        img_xray = drr_xray(rotation, translation_xray, parameterization=\"euler_angles\", convention=\"ZYX\").to(device)\n        img_mask = drr_mask(rotation, translation_mask, parameterization=\"euler_angles\", convention=\"ZYX\").to(device)\n        img_xray = img_xray.cpu()\n        img_xray = np.array(img_xray).squeeze()\n\n        img_mask = img_mask.cpu()\n        img_mask = np.array(img_mask).squeeze()\n\n#         patient_name_png = input_path_xray.split('/')[-2] + '.png'\n        patient_name_npy = input_path_xray.split('/')[-2] + '.npy'\n        \n        np.save(os.path.join(img_save_dir_xray, patient_name_npy), img_xray)\n        np.save(os.path.join(img_save_dir_mask, patient_name_npy), img_mask)\n\n#         plt.imsave(os.path.join(img_save_dir_png_xray, patient_name_png), img_xray, cmap='gray')\n#         plt.imsave(os.path.join(img_save_dir_png_mask, patient_name_png), img_mask, cmap='gray')\n    except:\n        print(\"error\")","metadata":{"execution":{"iopub.status.busy":"2024-01-15T07:26:55.319061Z","iopub.execute_input":"2024-01-15T07:26:55.319446Z","iopub.status.idle":"2024-01-15T07:27:58.264995Z","shell.execute_reply.started":"2024-01-15T07:26:55.31941Z","shell.execute_reply":"2024-01-15T07:27:58.263519Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# # Read in the volume and get the isocenter\n# volume, spacing = load_example_ct()\n# print(volume.shape)\n\n# plt.imshow(volume[:, :, 100], cmap='gray')\n# plt.show()\n# bx, by, bz = torch.tensor(volume.shape) * torch.tensor(spacing) / 2\n\n# # Initialize the DRR module for generating synthetic X-rays\n# device = \"cuda\" if torch.cuda.is_available() else \"cpu\"\n# drr = DRR(\n#     volume,      # The CT volume as a numpy array\n#     spacing,     # Voxel dimensions of the CT\n#     sdr=300.0,   # Source-to-detector radius (half of the source-to-detector distance)\n#     height=200,  # Height of the DRR (if width is not seperately provided, the generated image is square)\n#     delx=4.0,    # Pixel spacing (in mm)\n# ).to(device)\n\n# # Set the camera pose with rotation (yaw, pitch, roll) and translation (x, y, z)\n# rotation = torch.tensor([[torch.pi, 0.0, torch.pi / 2]], device=device)\n# translation = torch.tensor([[bx, by, bz]], device=device)\n\n# # 📸 Also note that DiffDRR can take many representations of SO(3) 📸\n# # For example, quaternions, rotation matrix, axis-angle, etc...\n# img = drr(rotation, translation, parameterization=\"euler_angles\", convention=\"ZYX\")\n# img = img.cpu()\n# img = np.array(img).squeeze()\n# plt.imsave(\"../tmp.png\", img, cmap='gray')\n# # plot_drr(img, ticks=False)\n# # plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# volume, spacing = load_nii_gz_file(\"../mask/abdominallymphnodes-00576/colon.nii.gz\")\n# print(\"Volume info: \", volume.shape, spacing)\n# # plt.imshow(volume[:, :, 100], cmap='gray')\n# # plt.show()\n# volume = np.transpose(volume, (2, 1, 0))\n# # アニメーションを作成\n# ani = create_animation(volume, \"TEST\")\n\n# # アニメーションを表示\n# HTML(ani.to_jshtml())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# volume, spacing = load_example_ct()\n# print(\"Volume info: \", volume.shape, spacing)\n# volume = np.transpose(volume, (2, 1, 0))\n# # アニメーションを作成\n# ani = create_animation(volume, \"TEST\")\n\n# # アニメーションを表示\n# HTML(ani.to_jshtml())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}