{"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":"# Implementation Example to fix PixelRepresentation == 1\n\nThe following will showcase two scenarios:\n1. Addressing the issue of visually unusual pixel arrays in DICOM files.\n2. Illustrate that the function does not affect DICOM files that appear visually normal.\n","metadata":{}},{"cell_type":"code","source":"# %pip install nnunetv2","metadata":{"execution":{"iopub.status.busy":"2023-08-03T13:07:10.642771Z","iopub.execute_input":"2023-08-03T13:07:10.643465Z","iopub.status.idle":"2023-08-03T13:07:10.66723Z","shell.execute_reply.started":"2023-08-03T13:07:10.643417Z","shell.execute_reply":"2023-08-03T13:07:10.665809Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd\nimport pydicom\nimport os\nimport gc\nfrom os.path import join, basename, dirname, exists\nfrom tqdm.notebook import tqdm\nfrom glob import glob\nimport numpy as np\nfrom matplotlib import pyplot as plt\nfrom scipy.ndimage import zoom\nfrom multiprocessing import Pool\n\nBASE_DIR = '/kaggle/input/rsna-2023-abdominal-trauma-detection'\n\ndef standardize_pixel_array(dcm: pydicom.dataset.FileDataset) -> np.ndarray:\n    # Correct DICOM pixel_array if PixelRepresentation == 1.\n    pixel_array = dcm.pixel_array\n    if dcm.PixelRepresentation == 1:\n        bit_shift = dcm.BitsAllocated - dcm.BitsStored\n        dtype = pixel_array.dtype \n        pixel_array = (pixel_array << bit_shift).astype(dtype) >>  bit_shift\n    return pixel_array # pydicom.pixel_data_handlers.util.apply_modality_lut(pixel_array, dcm)\n\n# train_dicom_tags = pd.read_parquet(join(BASE_DIR, 'train_dicom_tags.parquet'), engine='pyarrow')\n# %time all_paths = set(glob(join(BASE_DIR, 'train_images', '*', '*', '*.dcm')))\n# len(train_dicom_tags), len(all_paths), len(train_dicom_tags)-len(all_paths)","metadata":{"execution":{"iopub.status.busy":"2023-08-03T13:07:10.792786Z","iopub.execute_input":"2023-08-03T13:07:10.79346Z","iopub.status.idle":"2023-08-03T13:07:11.380007Z","shell.execute_reply.started":"2023-08-03T13:07:10.793418Z","shell.execute_reply":"2023-08-03T13:07:11.378794Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# %time len(glob(join(BASE_DIR, 'train_images', '*', '*')))","metadata":{"execution":{"iopub.status.busy":"2023-08-03T13:07:11.381839Z","iopub.execute_input":"2023-08-03T13:07:11.382192Z","iopub.status.idle":"2023-08-03T13:07:11.386625Z","shell.execute_reply.started":"2023-08-03T13:07:11.382162Z","shell.execute_reply":"2023-08-03T13:07:11.385447Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# %time lengths = [len(os.listdir(path)) for path in glob(join(BASE_DIR, 'train_images', '*', '*'))]\n# min(lengths), max(lengths)","metadata":{"execution":{"iopub.status.busy":"2023-08-03T13:07:11.387864Z","iopub.execute_input":"2023-08-03T13:07:11.388161Z","iopub.status.idle":"2023-08-03T13:07:11.397908Z","shell.execute_reply.started":"2023-08-03T13:07:11.388136Z","shell.execute_reply":"2023-08-03T13:07:11.396645Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# print(plt.hist(lengths))\n# plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-08-03T13:07:11.400458Z","iopub.execute_input":"2023-08-03T13:07:11.400786Z","iopub.status.idle":"2023-08-03T13:07:11.409728Z","shell.execute_reply.started":"2023-08-03T13:07:11.400758Z","shell.execute_reply":"2023-08-03T13:07:11.408548Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# diff_paths = set((BASE_DIR + '/' + train_dicom_tags.path).to_list()) - all_paths\n# len(diff_paths), 1510373-1500653","metadata":{"execution":{"iopub.status.busy":"2023-08-03T13:07:11.41126Z","iopub.execute_input":"2023-08-03T13:07:11.411621Z","iopub.status.idle":"2023-08-03T13:07:11.419725Z","shell.execute_reply.started":"2023-08-03T13:07:11.411593Z","shell.execute_reply":"2023-08-03T13:07:11.418832Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# missing_series = { basename(dirname(path)) for path in diff_paths }\n# len(missing_series), sum([exists(dirname(path)) for path in diff_paths])","metadata":{"execution":{"iopub.status.busy":"2023-08-03T13:07:11.420752Z","iopub.execute_input":"2023-08-03T13:07:11.421168Z","iopub.status.idle":"2023-08-03T13:07:11.432423Z","shell.execute_reply.started":"2023-08-03T13:07:11.421132Z","shell.execute_reply":"2023-08-03T13:07:11.431359Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# display(train_dicom_tags.Rows.value_counts())\n# display(train_dicom_tags.Columns.value_counts())","metadata":{"execution":{"iopub.status.busy":"2023-08-03T13:07:11.433495Z","iopub.execute_input":"2023-08-03T13:07:11.433897Z","iopub.status.idle":"2023-08-03T13:07:11.443926Z","shell.execute_reply.started":"2023-08-03T13:07:11.433861Z","shell.execute_reply":"2023-08-03T13:07:11.442907Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# print(train_dicom_tags[train_dicom_tags.Columns != train_dicom_tags.Rows].Rows.unique())\n# train_dicom_tags[train_dicom_tags.Columns != train_dicom_tags.Rows].Columns.value_counts().sort_index()","metadata":{"execution":{"iopub.status.busy":"2023-08-03T13:07:11.444874Z","iopub.execute_input":"2023-08-03T13:07:11.445217Z","iopub.status.idle":"2023-08-03T13:07:11.455873Z","shell.execute_reply.started":"2023-08-03T13:07:11.44518Z","shell.execute_reply":"2023-08-03T13:07:11.454856Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# train_dicom_tags.columns #PatientID,InstanceNumber,SeriesNumber,StudyInstanceUID\n# display(train_dicom_tags[train_dicom_tags.path=='train_images/59261/5980/181.dcm'].T) # [['PatientID','SeriesInstanceUID','InstanceNumber']])\n# len(train_dicom_tags[train_dicom_tags.SeriesInstanceUID.str.endswith('.5980')])","metadata":{"execution":{"iopub.status.busy":"2023-08-03T13:07:11.457256Z","iopub.execute_input":"2023-08-03T13:07:11.457724Z","iopub.status.idle":"2023-08-03T13:07:11.466666Z","shell.execute_reply.started":"2023-08-03T13:07:11.457685Z","shell.execute_reply":"2023-08-03T13:07:11.465733Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_volume_from_dcm_folder(directory):\n    return np.stack([\n        standardize_pixel_array(pydicom.read_file(path))\n        for path in sorted(glob(join(directory, '*.dcm')))\n    ])\n\n# %time vol_arr = get_volume_from_dcm_folder(join(BASE_DIR, 'train_images/10004/21057'))\n# vol_arr.shape","metadata":{"execution":{"iopub.status.busy":"2023-08-03T13:07:11.471423Z","iopub.execute_input":"2023-08-03T13:07:11.471839Z","iopub.status.idle":"2023-08-03T13:07:11.48671Z","shell.execute_reply.started":"2023-08-03T13:07:11.471764Z","shell.execute_reply":"2023-08-03T13:07:11.485615Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# d, h, w = vol_arr.shape\n# %time zoom(vol_arr, (256/d,512/h,512/w), order=3).shape","metadata":{"execution":{"iopub.status.busy":"2023-08-03T13:07:11.488212Z","iopub.execute_input":"2023-08-03T13:07:11.48876Z","iopub.status.idle":"2023-08-03T13:07:11.498316Z","shell.execute_reply.started":"2023-08-03T13:07:11.488726Z","shell.execute_reply":"2023-08-03T13:07:11.497147Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# !rm *.npy\n\n# for dir_path in tqdm(glob(join(BASE_DIR, 'train_images', '*', '*'))):\n#     vol_arr = get_volume_from_dcm_folder(dir_path)\n#     d, h, w = vol_arr.shape\n#     vol_arr = zoom(vol_arr, (128/d,128/h,128/w), order=3)\n#     np.save(f\"{basename(dir_path)}.npy\", vol_arr)","metadata":{"execution":{"iopub.status.busy":"2023-08-03T13:07:11.499874Z","iopub.execute_input":"2023-08-03T13:07:11.500219Z","iopub.status.idle":"2023-08-03T13:07:11.510948Z","shell.execute_reply.started":"2023-08-03T13:07:11.500182Z","shell.execute_reply":"2023-08-03T13:07:11.509962Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def f(dir_path):\n    vol_arr = get_volume_from_dcm_folder(dir_path)\n    d, h, w = vol_arr.shape\n    vol_arr = zoom(vol_arr, (128/d,128/h,128/w), order=3)\n    filename = f\"{basename(dir_path)}.npy\"\n    np.save(filename, vol_arr)\n    del vol_arr\n    # pbar.update()\n    gc.collect()\n    print(f\"Done with {filename}\\n\", end=\"\")\n    \n\ndir_paths = glob(join(BASE_DIR, 'train_images', '*', '*'))\n# with tqdm(total=len(dir_paths)) as pbar:\nif __name__ == '__main__':\n    with Pool(5) as p:\n        print(p.map(f, dir_paths))","metadata":{"execution":{"iopub.status.busy":"2023-08-03T13:07:11.511933Z","iopub.execute_input":"2023-08-03T13:07:11.512347Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Example Case of Unusual DICOM Pixel Array\n\nExample Case of Unusual DICOM Pixel Array\npatient_id = 5425\n\nExample image is found in index 917637 of pixel_representation pandas DataFrame.\n\nYou will find that the Original image appears very unusual, but fixed after implementing the function.","metadata":{}},{"cell_type":"code","source":"# # Grab the path of the image in index 917637\n# sample_image = train_dicom_tags.loc[917637, 'path']\n\n# # Open the DICOM file using pydicom.\n# dcm = pydicom.read_file(join(BASE_DIR, sample_image))\n\n# original_pixel_array = dcm.pixel_array\n# fixed_pixel_array = standardize_pixel_array(dcm)\n\n# # To demonstrate that Pixel Repsentation is indeed 1\n# print(f'{dcm.SOPInstanceUID}: Pixel Representation = {dcm.PixelRepresentation}')\n# # To demonstrate whether the pixel array has been changed after passing through the function.\n# print(f'Two images are equal: {np.array_equal(original_pixel_array,fixed_pixel_array)}')\n\n# # plot both images\n# fig, (ax1, ax2) = plt.subplots(1, 2)\n# ax1.imshow(original_pixel_array, cmap='gray')\n# ax1.set_title('Original')\n# ax1.axis('off')\n# ax2.imshow(fixed_pixel_array, cmap='gray')\n# ax2.set_title('Fixed')\n# ax2.axis('off')\n# fig.tight_layout()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Example Case of Normal DICOM Pixel Array\npatient_id = 19\n\nExample image is found in index 1477490 of pixel_representation pandas DataFrame.\n\nYou will find that even after the DICOM pixel array has undergone processing through the function, its content array unaltered if the DICOM image appears normal.\n","metadata":{}},{"cell_type":"code","source":"# # Grab the path of the image in index 1477490\n# sample_image = train_dicom_tags.loc[1477490]['path']\n\n# # Open the DICOM file using pydicom.\n# dcm = pydicom.read_file(join(BASE_DIR, sample_image))\n\n# original_pixel_array = dcm.pixel_array\n# fixed_pixel_array = standardize_pixel_array(dcm)\n\n# # To demonstrate that Pixel Repsentation is indeed 1\n# print(f'{dcm.SOPInstanceUID}: Pixel Representation = {dcm.PixelRepresentation}')\n# # To demonstrate whether the pixel array has been changed after passing through the function.\n# print(f'Two images are equal: {np.array_equal(original_pixel_array,fixed_pixel_array)}')\n\n# # plot both images\n# fig, (ax1, ax2) = plt.subplots(1, 2)\n# ax1.imshow(original_pixel_array, cmap='gray')\n# ax1.set_title('Original')\n# ax1.axis('off')\n# ax2.imshow(fixed_pixel_array, cmap='gray')\n# ax2.set_title('Fixed')\n# ax2.axis('off')\n# fig.tight_layout()","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}