{"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":"# Resampling pipeline for RSNA 2023 Abdominal Trauma Detection\n\nThis notebook loads and resamples the original images data to the same spacing (2mm, 2mm, 2mm) and size(320, 256, 256) and builds datasets and dataloaders with PyTorch.\n\nThanks to the functions to load .dcm files and construct 3D images provided from [https://www.kaggle.com/code/parhammostame/construct-3d-arrays-from-dcm-nii-3-view-angles](http://)\n\nUsing dataloader num_workers = 4, the preprocessing of test set can be done in about 1h 10 mins","metadata":{}},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\n# for dirname, _, filenames in os.walk('/kaggle/input'):\n#     for filename in filenames:\n#         print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:01:51.841578Z","iopub.execute_input":"2023-08-20T22:01:51.841943Z","iopub.status.idle":"2023-08-20T22:01:51.848357Z","shell.execute_reply.started":"2023-08-20T22:01:51.841908Z","shell.execute_reply":"2023-08-20T22:01:51.847382Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nimport os\nimport random\nimport re\n\nfrom tqdm import tqdm\n\nimport pydicom as dicom\nimport nibabel as nib\nimport SimpleITK as sitk\n\nimport torch\nimport torch.nn as nn","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:01:51.853003Z","iopub.execute_input":"2023-08-20T22:01:51.853459Z","iopub.status.idle":"2023-08-20T22:01:51.867111Z","shell.execute_reply.started":"2023-08-20T22:01:51.853433Z","shell.execute_reply":"2023-08-20T22:01:51.866089Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"SEED = 344\ndevice = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\ndef seed_everything(seed):\n    random.seed(seed)\n    os.environ['PYTHONHASHSEED'] = str(seed)\n    np.random.seed(seed)\n    torch.manual_seed(seed)\n    torch.cuda.manual_seed(seed)\n    torch.backends.cudnn.deterministic = True # Fix the network according to random seed\n    print('Finish seeding with seed {}'.format(seed))\n    \nseed_everything(SEED)\nprint('Training on device {}'.format(device))\n\nimport multiprocessing\nprint(multiprocessing.cpu_count())","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:01:51.869251Z","iopub.execute_input":"2023-08-20T22:01:51.869952Z","iopub.status.idle":"2023-08-20T22:01:51.884347Z","shell.execute_reply.started":"2023-08-20T22:01:51.869918Z","shell.execute_reply":"2023-08-20T22:01:51.883314Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dicom_tag_columns = [\n    'Columns',\n    'ImageOrientationPatient',\n    'ImagePositionPatient',\n    'InstanceNumber',\n    'PatientID',\n    'PatientPosition',\n    'PixelSpacing',\n    'RescaleIntercept',\n    'RescaleSlope',\n    'Rows',\n    'SeriesNumber',\n    'SliceThickness',\n    'path',\n    'WindowCenter',\n    'WindowWidth'\n]\n\ntrain_dicom_tags = pd.read_parquet('/kaggle/input/rsna-2023-abdominal-trauma-detection/train_dicom_tags.parquet', columns=dicom_tag_columns)\ntest_dicom_tags = pd.read_parquet('/kaggle/input/rsna-2023-abdominal-trauma-detection/test_dicom_tags.parquet', columns=dicom_tag_columns)\n\ntrain_series_meta = pd.read_csv('/kaggle/input/rsna-2023-abdominal-trauma-detection/train_series_meta.csv')\ntest_series_meta = pd.read_csv('/kaggle/input/rsna-2023-abdominal-trauma-detection/test_series_meta.csv')\n\ntest_dicom_tags","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:01:51.885626Z","iopub.execute_input":"2023-08-20T22:01:51.886662Z","iopub.status.idle":"2023-08-20T22:01:54.091422Z","shell.execute_reply.started":"2023-08-20T22:01:51.88654Z","shell.execute_reply":"2023-08-20T22:01:54.090448Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#train_dicom_tags[train_dicom_tags.RescaleIntercept == 0]\ntrain_dicom_tags[train_dicom_tags.RescaleIntercept == 0]","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:01:54.094322Z","iopub.execute_input":"2023-08-20T22:01:54.094711Z","iopub.status.idle":"2023-08-20T22:01:54.205556Z","shell.execute_reply.started":"2023-08-20T22:01:54.094676Z","shell.execute_reply":"2023-08-20T22:01:54.204354Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dfilter = (train_dicom_tags['SliceThickness'] == 3.0)\ndfilter = train_dicom_tags[dfilter].PatientID.drop_duplicates().values.astype('int64')\ntrain_series_meta[train_series_meta['patient_id'].isin(dfilter)]","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:01:54.207552Z","iopub.execute_input":"2023-08-20T22:01:54.207969Z","iopub.status.idle":"2023-08-20T22:01:54.274834Z","shell.execute_reply.started":"2023-08-20T22:01:54.207932Z","shell.execute_reply":"2023-08-20T22:01:54.273817Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def raw_path_gen(patient_id, series_id, train=True):\n    path = '/kaggle/input/rsna-2023-abdominal-trauma-detection/'\n    if(train):\n        path += 'train_images/'\n    else:\n        path += 'test_images/'\n    \n    path += str(patient_id) + '/' + str(series_id)\n    \n    return path\n\ndef create_3D_scans(folder, downsample_rate=1): \n    filenames = os.listdir(folder)\n    filenames = [int(filename.split('.')[0]) for filename in filenames]\n    filenames = sorted(filenames)\n    filenames = [str(filename) + '.dcm' for filename in filenames]\n        \n    volume = []\n    #for filename in tqdm(filenames[::downsample_rate], position=0): \n    for filename in filenames[::downsample_rate]: \n        filepath = os.path.join(folder, filename)\n        ds = dicom.dcmread(filepath)\n        image = ds.pixel_array\n\n        \n        if ds.PixelRepresentation == 1:\n            bit_shift = ds.BitsAllocated - ds.BitsStored\n            dtype = image.dtype \n            image = (image << bit_shift).astype(dtype) >>  bit_shift\n        \n        # find rescale params\n        if (\"RescaleIntercept\" in ds) and (\"RescaleSlope\" in ds):\n            intercept = float(ds.RescaleIntercept)\n            slope = float(ds.RescaleSlope)\n    \n        # find clipping params\n        center = int(ds.WindowCenter)\n        width = int(ds.WindowWidth)\n        low = center - width / 2\n        high = center + width / 2    \n        \n        \n        image = (image * slope) + intercept\n        image = np.clip(image, low, high)\n\n        image = (image / np.max(image) * 255).astype(np.int16)\n        image = image[::downsample_rate, ::downsample_rate]\n        volume.append( image )\n    \n    volume = np.stack(volume, axis=0)\n    return volume","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:01:54.276606Z","iopub.execute_input":"2023-08-20T22:01:54.277325Z","iopub.status.idle":"2023-08-20T22:01:54.294186Z","shell.execute_reply.started":"2023-08-20T22:01:54.277284Z","shell.execute_reply":"2023-08-20T22:01:54.293135Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def test_3D_scans(folder, downsample_rate=1): \n    filenames = os.listdir(folder)\n    filenames = [int(filename.split('.')[0]) for filename in filenames]\n    filenames = sorted(filenames)\n    filenames = [str(filename) + '.dcm' for filename in filenames]\n        \n    volume = []\n    #for filename in tqdm(filenames[::downsample_rate], position=0): \n    for filename in filenames[::downsample_rate]: \n        filepath = os.path.join(folder, filename)\n        ds = dicom.dcmread(filepath)\n        image = ds.pixel_array\n        \n        if ds.PixelRepresentation == 1:\n            print(filepath)\n            bit_shift = ds.BitsAllocated - ds.BitsStored\n            dtype = image.dtype \n            image = (image << bit_shift).astype(dtype) >>  bit_shift\n            return\n        break","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:01:54.295673Z","iopub.execute_input":"2023-08-20T22:01:54.296282Z","iopub.status.idle":"2023-08-20T22:01:54.316715Z","shell.execute_reply.started":"2023-08-20T22:01:54.296245Z","shell.execute_reply":"2023-08-20T22:01:54.315874Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# for i in tqdm(range(0, len(train_series_meta)), position=0):\n#     patient_id, series_id = train_series_meta.loc[i, [\"patient_id\", \"series_id\"]].astype('int')\n#     filepath = raw_path_gen(patient_id, series_id)\n#     test_3D_scans(filepath)","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:01:54.318127Z","iopub.execute_input":"2023-08-20T22:01:54.318721Z","iopub.status.idle":"2023-08-20T22:01:54.333346Z","shell.execute_reply.started":"2023-08-20T22:01:54.318683Z","shell.execute_reply":"2023-08-20T22:01:54.33228Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_image_with_seg(volume, volume_seg=[], orientation='Coronal', num_subplots=20):\n    # simply copy\n    if len(volume_seg) == 0:\n        plot_mask = 0\n    else:\n        plot_mask = 1\n        \n    if orientation == 'Coronal':\n        slices = np.linspace(0, volume.shape[2]-1, num_subplots).astype(np.int16)\n        volume = volume.transpose([1, 0, 2])\n        if plot_mask:\n            volume_seg = volume_seg.transpose([1, 0, 2])\n        \n    elif orientation == 'Sagittal':\n        slices = np.linspace(0, volume.shape[2]-1, num_subplots).astype(np.int16)\n        volume = volume.transpose([2, 0, 1])\n        if plot_mask:\n            volume_seg = volume_seg.transpose([2, 0, 1])\n\n    elif orientation == 'Axial':\n        slices = np.linspace(0, volume.shape[0]-1, num_subplots).astype(np.int16)\n           \n    rows = np.max( [np.floor(np.sqrt(num_subplots)).astype(int) - 2, 1])\n    cols = np.ceil(num_subplots/rows).astype(int)\n    \n    fig, ax = plt.subplots(rows, cols, figsize=(cols * 2, rows * 4))\n    fig.tight_layout(h_pad=0.01, w_pad=0)\n    \n    ax = ax.ravel()\n    for this_ax in ax:\n        this_ax.axis('off')\n\n    for counter, this_slice in enumerate( slices ):\n        plt.sca(ax[counter])\n        \n        image = volume[this_slice, :, :]\n        plt.imshow(image, cmap='gray')\n        \n        if plot_mask:\n            mask = np.where(volume_seg[this_slice, :, :], volume_seg[this_slice, :, :], np.nan)\n            plt.imshow(mask, cmap='Set1', alpha=0.5)  ","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:01:54.335049Z","iopub.execute_input":"2023-08-20T22:01:54.336104Z","iopub.status.idle":"2023-08-20T22:01:54.35283Z","shell.execute_reply.started":"2023-08-20T22:01:54.335918Z","shell.execute_reply":"2023-08-20T22:01:54.351901Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# filepath = raw_path_gen(14429, 57624)\n# volume = create_3D_scans(filepath)\n# print(f'3D Image file shape: {volume.shape}')","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:01:54.357975Z","iopub.execute_input":"2023-08-20T22:01:54.358298Z","iopub.status.idle":"2023-08-20T22:01:54.36887Z","shell.execute_reply.started":"2023-08-20T22:01:54.358273Z","shell.execute_reply":"2023-08-20T22:01:54.367885Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# plot_image_with_seg(volume, orientation='Axial', num_subplots=5)\n# plot_image_with_seg(volume, orientation='Sagittal', num_subplots=5)","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:01:54.370421Z","iopub.execute_input":"2023-08-20T22:01:54.370877Z","iopub.status.idle":"2023-08-20T22:01:54.380861Z","shell.execute_reply.started":"2023-08-20T22:01:54.370842Z","shell.execute_reply":"2023-08-20T22:01:54.37989Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#load series and resample\nimport re\n\nfake_img = sitk.ReadImage('/kaggle/input/rsna-2023-abdominal-trauma-detection/test_images/48843/62825/30.dcm')\n\ndef load_resample(patient_id, series_id, train = True,\n                  target_spacing = (2.0, 2.0, 2.0),\n                  target_size = (256, 256, 320)):\n    filepath = raw_path_gen(patient_id, series_id, train=train)\n    if not os.path.exists(filepath):\n        return fake_img\n    \n    if (train):\n        dicom_tags = train_dicom_tags\n    else:\n        dicom_tags = test_dicom_tags\n    \n    filenames = os.listdir(filepath)\n    \n    if(len(filenames) < 2):\n        return sitk.GetImageFromArray(create_3D_scans(filepath))\n    \n    dicom_name1 = filepath + '/' + filenames[0]\n    dicom_name1 = dicom_name1.split('/')[-4:]\n    dicom_name1 = '/'.join(dicom_name1)\n    dicom_tags1 = dicom_tags[dicom_tags.path == dicom_name1]\n    \n    dicom_name2 = filepath + '/' + filenames[1]\n    dicom_name2 = dicom_name2.split('/')[-4:]\n    dicom_name2 = '/'.join(dicom_name2)\n    dicom_tags2 = dicom_tags[dicom_tags.path == dicom_name2]\n    \n    dicom_po1 = dicom_tags1.ImagePositionPatient.values[0]\n    dicom_po2 = dicom_tags2.ImagePositionPatient.values[0]\n    dicom_po1 = float(re.split(', |]', dicom_po1)[-2])\n    dicom_po2 = float(re.split(', |]', dicom_po2)[-2])\n    \n    dicom_dz = abs((dicom_po1 - dicom_po2) / (dicom_tags1.InstanceNumber.values[0] - dicom_tags2.InstanceNumber.values[0]))\n\n    dicom_spacing_x = float(re.split(', |]|\\[', dicom_tags1.PixelSpacing.values[0])[-2])\n    dicom_spacing_y = float(re.split(', |]|\\[', dicom_tags1.PixelSpacing.values[0])[1])\n\n    if (dicom_dz <= 1.0):\n        downsample_rate = 2\n    else:\n        downsample_rate = 1\n    \n    dicom_spacing_x = dicom_spacing_x * downsample_rate\n    dicom_spacing_y = dicom_spacing_y * downsample_rate\n    dicom_dz = dicom_dz * downsample_rate\n    \n    original_image = create_3D_scans(filepath, downsample_rate)\n    original_image = sitk.GetImageFromArray(original_image)\n    original_image.SetSpacing((dicom_spacing_x,\n                        dicom_spacing_y,\n                        dicom_dz))\n    \n\n    original_size = original_image.GetSize()\n    original_spacing = original_image.GetSpacing()\n    \n    r_x = dicom_spacing_x / target_spacing[0]\n    r_y = dicom_spacing_y / target_spacing[1]\n    r_z = dicom_dz / target_spacing[2]\n    \n    o_x = -(target_size[0] - original_size[0] * r_x)\n    o_y = -(target_size[0] - original_size[1] * r_y)\n    o_z = -(target_size[2] - original_size[2] * r_z)\n    #print((target_size[0] - original_size[0] * r_x) / 2)\n    \n    \n    target_origin = (o_x, o_y, o_z)\n    \n    \n    #print(original_spacing)\n    #print(original_size)\n    \n    resampler = sitk.ResampleImageFilter()\n    resampler.SetInterpolator(sitk.sitkLinear)\n    resampler.SetOutputSpacing(target_spacing)\n    resampler.SetSize(target_size)\n    resampler.SetOutputOrigin(target_origin)\n    resampler.SetDefaultPixelValue(-150)\n    resampled_image = resampler.Execute(original_image)\n    #print(resampled_image.GetSize())\n    \n    \n    return resampled_image\n\n    #49954, 41479\n    #10026, 42932\n    \n    \n# img = load_resample(10004, 21057)\n# sitk.WriteImage(img, './test.nii.gz')\n# image_a = sitk.GetArrayFromImage(img)\n# #image_a = np.flip(image_a, 0)\n# #image_a = np.clip(image_a, -150, 250)\n\n# print(image_a.shape)\n# plot_image_with_seg(image_a, orientation='Axial', num_subplots=5)\n# plot_image_with_seg(image_a, orientation='Sagittal', num_subplots=5)","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:01:54.383973Z","iopub.execute_input":"2023-08-20T22:01:54.384347Z","iopub.status.idle":"2023-08-20T22:01:54.498012Z","shell.execute_reply.started":"2023-08-20T22:01:54.384322Z","shell.execute_reply":"2023-08-20T22:01:54.496972Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from torch.utils.data import Dataset\nfrom torch.utils.data import DataLoader\n\nclass CTDataset(Dataset):\n    def __init__(self, train=True, device='cpu'):\n        self.train = train\n        self.device = device\n        if(train):\n            self.series_meta = train_series_meta\n        else:\n            self.series_meta = test_series_meta\n        \n    def __len__(self):\n        #return 1100\n        return len(self.series_meta)\n    \n    def __getitem__(self, idx):\n        patient_id, series_id = self.series_meta.loc[idx, [\"patient_id\", \"series_id\"]].astype('int')\n        img_a = sitk.GetArrayFromImage(load_resample(patient_id, series_id, train=self.train))\n        img_t = torch.from_numpy(img_a).to(self.device)\n        return img_t","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:01:54.500252Z","iopub.execute_input":"2023-08-20T22:01:54.500947Z","iopub.status.idle":"2023-08-20T22:01:54.509463Z","shell.execute_reply.started":"2023-08-20T22:01:54.500909Z","shell.execute_reply":"2023-08-20T22:01:54.508331Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import multiprocessing\nnum_cpus = multiprocessing.cpu_count()\n\ntrain_ds = CTDataset(train=True)\ntrain_dl = DataLoader(train_ds, batch_size=4, shuffle=True, num_workers=num_cpus)\n\ntest_ds = CTDataset(train=False)\ntest_dl = DataLoader(test_ds, batch_size=4, shuffle=True, num_workers=num_cpus)","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:01:54.510974Z","iopub.execute_input":"2023-08-20T22:01:54.511349Z","iopub.status.idle":"2023-08-20T22:01:54.526269Z","shell.execute_reply.started":"2023-08-20T22:01:54.511294Z","shell.execute_reply":"2023-08-20T22:01:54.525282Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for imgs in tqdm(test_dl):\n    print(imgs.shape)","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:01:54.527939Z","iopub.execute_input":"2023-08-20T22:01:54.52839Z","iopub.status.idle":"2023-08-20T22:01:54.736902Z","shell.execute_reply.started":"2023-08-20T22:01:54.528357Z","shell.execute_reply":"2023-08-20T22:01:54.735693Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# img_a = sitk.GetArrayFromImage(load_resample(49954, 41479, train=True))\n# plot_image_with_seg(img_a, orientation='Axial', num_subplots=5)\n# plot_image_with_seg(img_a, orientation='Sagittal', num_subplots=5)","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:01:54.739063Z","iopub.execute_input":"2023-08-20T22:01:54.739459Z","iopub.status.idle":"2023-08-20T22:01:54.745398Z","shell.execute_reply.started":"2023-08-20T22:01:54.739417Z","shell.execute_reply":"2023-08-20T22:01:54.743332Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!cp /kaggle/input/rsna-2023-abdominal-trauma-detection/sample_submission.csv submission.csv","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:01:54.746769Z","iopub.execute_input":"2023-08-20T22:01:54.747555Z","iopub.status.idle":"2023-08-20T22:01:55.769926Z","shell.execute_reply.started":"2023-08-20T22:01:54.747518Z","shell.execute_reply":"2023-08-20T22:01:55.768567Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}