{"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":"code","source":"import os\nimport sys \nimport json\nimport glob\nimport random\nimport collections\nimport time\nimport re\n\nimport numpy as np\nimport pandas as pd\nimport pydicom\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\nimport cv2\nimport PIL.Image\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nfrom tqdm import tqdm","metadata":{"papermill":{"duration":1.048295,"end_time":"2021-07-14T20:26:46.309722","exception":false,"start_time":"2021-07-14T20:26:45.261427","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-08-16T21:55:56.399641Z","iopub.execute_input":"2021-08-16T21:55:56.399987Z","iopub.status.idle":"2021-08-16T21:55:56.406672Z","shell.execute_reply.started":"2021-08-16T21:55:56.399958Z","shell.execute_reply":"2021-08-16T21:55:56.405680Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data_directory = '../input/rsna-miccai-brain-tumor-radiogenomic-classification'","metadata":{"lines_to_end_of_cell_marker":2,"lines_to_next_cell":2,"papermill":{"duration":0.05565,"end_time":"2021-07-14T20:26:46.486521","exception":false,"start_time":"2021-07-14T20:26:46.430871","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2021-08-16T21:55:56.418217Z","iopub.execute_input":"2021-08-16T21:55:56.418711Z","iopub.status.idle":"2021-08-16T21:55:56.422073Z","shell.execute_reply.started":"2021-08-16T21:55:56.418673Z","shell.execute_reply":"2021-08-16T21:55:56.421346Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Set up your config","metadata":{}},{"cell_type":"code","source":"DEBUG = False\nmri_types = ['FLAIR','T1w','T1wCE','T2w']\nSIZE = 256\nNUM_IMAGES = 64","metadata":{"execution":{"iopub.status.busy":"2021-08-16T21:55:56.450836Z","iopub.execute_input":"2021-08-16T21:55:56.451350Z","iopub.status.idle":"2021-08-16T21:55:56.454987Z","shell.execute_reply.started":"2021-08-16T21:55:56.451310Z","shell.execute_reply":"2021-08-16T21:55:56.454294Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Let's read in one dicom file to show you what I'm talking about.","metadata":{}},{"cell_type":"markdown","source":"# Scenario 1: Read it WITHOUT applying voi_lut","metadata":{}},{"cell_type":"code","source":"files = sorted(glob.glob(f\"{data_directory}/train/00688/FLAIR/*.dcm\"), \n               key=lambda var:[int(x) if x.isdigit() else x for x in re.findall(r'[^0-9]|[0-9]+', var)])","metadata":{"execution":{"iopub.status.busy":"2021-08-16T21:55:56.471874Z","iopub.execute_input":"2021-08-16T21:55:56.472453Z","iopub.status.idle":"2021-08-16T21:55:56.490403Z","shell.execute_reply.started":"2021-08-16T21:55:56.472410Z","shell.execute_reply":"2021-08-16T21:55:56.489699Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"slices = []\nfor filepath in tqdm(files):\n    dicom = pydicom.read_file(filepath)\n    data = dicom.pixel_array\n    slices.append(data)\n    \n# Also, get the plane of the DICOM - axial, sagittal, or coronal.\nx1, y1, _, x2, y2, _ = [round(j) for j in dicom.ImageOrientationPatient]\ncords = [x1, y1, x2, y2]\n\nif cords == [1, 0, 0, 0]:\n    plane = 'Coronal'\nelif cords == [1, 0, 0, 1]:\n    plane = 'Axial'\nelif cords == [0, 1, 0, 0]:\n    plane = 'Sagittal'\nelse:\n    plane = 'Unknown'\nprint(plane)","metadata":{"execution":{"iopub.status.busy":"2021-08-16T21:55:56.491787Z","iopub.execute_input":"2021-08-16T21:55:56.492319Z","iopub.status.idle":"2021-08-16T21:55:56.876752Z","shell.execute_reply.started":"2021-08-16T21:55:56.492276Z","shell.execute_reply":"2021-08-16T21:55:56.876077Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"img3d = np.stack(slices)","metadata":{"execution":{"iopub.status.busy":"2021-08-16T21:55:56.877930Z","iopub.execute_input":"2021-08-16T21:55:56.878303Z","iopub.status.idle":"2021-08-16T21:55:56.914952Z","shell.execute_reply.started":"2021-08-16T21:55:56.878267Z","shell.execute_reply":"2021-08-16T21:55:56.914252Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Normalize it so we can plot it\nif img3d.sum() != 0:\n    img3d = img3d - np.min(img3d)\n    img3d = img3d / np.max(img3d)\n    img3d = (img3d * 255).astype(np.uint8)","metadata":{"execution":{"iopub.status.busy":"2021-08-16T21:55:56.916142Z","iopub.execute_input":"2021-08-16T21:55:56.916550Z","iopub.status.idle":"2021-08-16T21:55:57.426909Z","shell.execute_reply.started":"2021-08-16T21:55:56.916499Z","shell.execute_reply":"2021-08-16T21:55:57.425978Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(img3d.shape)","metadata":{"execution":{"iopub.status.busy":"2021-08-16T21:55:57.428193Z","iopub.execute_input":"2021-08-16T21:55:57.428606Z","iopub.status.idle":"2021-08-16T21:55:57.433715Z","shell.execute_reply.started":"2021-08-16T21:55:57.428566Z","shell.execute_reply":"2021-08-16T21:55:57.432796Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# The coronal view of the scan\nplt.imshow(img3d[len(img3d)//2, :, :])","metadata":{"execution":{"iopub.status.busy":"2021-08-16T21:55:57.435308Z","iopub.execute_input":"2021-08-16T21:55:57.435766Z","iopub.status.idle":"2021-08-16T21:55:57.763313Z","shell.execute_reply.started":"2021-08-16T21:55:57.435727Z","shell.execute_reply":"2021-08-16T21:55:57.762433Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# The axial view of the scan\nplt.imshow(img3d[:, img3d.shape[1]//2, :])","metadata":{"execution":{"iopub.status.busy":"2021-08-16T21:56:19.691340Z","iopub.execute_input":"2021-08-16T21:56:19.691693Z","iopub.status.idle":"2021-08-16T21:56:19.843093Z","shell.execute_reply.started":"2021-08-16T21:56:19.691661Z","shell.execute_reply":"2021-08-16T21:56:19.842259Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Quickly plot max values for each slice on the axial view","metadata":{}},{"cell_type":"code","source":"maxs = []\nq25s = []\nfor slice_idx in tqdm(range(len(img3d))):\n    maxs.append(np.max(img3d[slice_idx, img3d.shape[1]//2, :]))\n    q25s.append(np.quantile(img3d[slice_idx, img3d.shape[1]//2, :], 0.25))\n    \nax = sns.lineplot(x = range(len(img3d)), y = maxs, color='orange', label='Max values')\nsns.lineplot(x = range(len(img3d)), y=q25s, color='blue', label='Q25 values')","metadata":{"execution":{"iopub.status.busy":"2021-08-16T21:57:03.238083Z","iopub.execute_input":"2021-08-16T21:57:03.238411Z","iopub.status.idle":"2021-08-16T21:57:03.503740Z","shell.execute_reply.started":"2021-08-16T21:57:03.238383Z","shell.execute_reply":"2021-08-16T21:57:03.502804Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Scenario 2: Repeat the same code, but applying VOI LUT.","metadata":{}},{"cell_type":"code","source":"slices = []\nfor filepath in tqdm(files):\n    dicom = pydicom.read_file(filepath)\n    data = apply_voi_lut(dicom.pixel_array, dicom)\n    slices.append(data)\n    \n# Also, get the plane of the DICOM - axial, sagittal, or coronal.\nx1, y1, _, x2, y2, _ = [round(j) for j in dicom.ImageOrientationPatient]\ncords = [x1, y1, x2, y2]\n\nif cords == [1, 0, 0, 0]:\n    plane = 'Coronal'\nelif cords == [1, 0, 0, 1]:\n    plane = 'Axial'\nelif cords == [0, 1, 0, 0]:\n    plane = 'Sagittal'\nelse:\n    plane = 'Unknown'\nprint(plane)","metadata":{"execution":{"iopub.status.busy":"2021-08-16T21:57:09.582332Z","iopub.execute_input":"2021-08-16T21:57:09.582656Z","iopub.status.idle":"2021-08-16T21:57:10.540465Z","shell.execute_reply.started":"2021-08-16T21:57:09.582628Z","shell.execute_reply":"2021-08-16T21:57:10.539558Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"img3d = np.stack(slices)\n# Normalize it so we can plot it\nif img3d.sum() != 0:\n    img3d = img3d - np.min(img3d)\n    img3d = img3d / np.max(img3d)\n    img3d = (img3d * 255).astype(np.uint8)","metadata":{"execution":{"iopub.status.busy":"2021-08-16T21:57:11.940421Z","iopub.execute_input":"2021-08-16T21:57:11.940755Z","iopub.status.idle":"2021-08-16T21:57:12.606194Z","shell.execute_reply.started":"2021-08-16T21:57:11.940728Z","shell.execute_reply":"2021-08-16T21:57:12.605370Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(img3d.shape)","metadata":{"execution":{"iopub.status.busy":"2021-08-16T21:57:12.607484Z","iopub.execute_input":"2021-08-16T21:57:12.607762Z","iopub.status.idle":"2021-08-16T21:57:12.612507Z","shell.execute_reply.started":"2021-08-16T21:57:12.607736Z","shell.execute_reply":"2021-08-16T21:57:12.611523Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# The coronal view of the scan\nplt.imshow(img3d[len(img3d)//2, :, :])","metadata":{"execution":{"iopub.status.busy":"2021-08-16T21:57:12.613971Z","iopub.execute_input":"2021-08-16T21:57:12.614242Z","iopub.status.idle":"2021-08-16T21:57:12.794731Z","shell.execute_reply.started":"2021-08-16T21:57:12.614216Z","shell.execute_reply":"2021-08-16T21:57:12.793783Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# The axial view of the scan\nplt.imshow(img3d[:, img3d.shape[1]//2, :])","metadata":{"execution":{"iopub.status.busy":"2021-08-16T21:57:27.674437Z","iopub.execute_input":"2021-08-16T21:57:27.674810Z","iopub.status.idle":"2021-08-16T21:57:27.820816Z","shell.execute_reply.started":"2021-08-16T21:57:27.674781Z","shell.execute_reply":"2021-08-16T21:57:27.819962Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Quickly plot max values for each slice","metadata":{}},{"cell_type":"code","source":"maxs = []\nq25s = []\nfor slice_idx in tqdm(range(len(img3d))):\n    maxs.append(np.max(img3d[slice_idx, img3d.shape[1]//2, :]))\n    q25s.append(np.quantile(img3d[slice_idx, img3d.shape[1]//2, :], 0.25))\n    \nax = sns.lineplot(x = range(len(img3d)), y = maxs, color='orange', label='Max values')\nsns.lineplot(x = range(len(img3d)), y=q25s, color='blue', label='Q25 values')","metadata":{"execution":{"iopub.status.busy":"2021-08-16T21:57:37.428615Z","iopub.execute_input":"2021-08-16T21:57:37.429000Z","iopub.status.idle":"2021-08-16T21:57:37.680420Z","shell.execute_reply.started":"2021-08-16T21:57:37.428969Z","shell.execute_reply":"2021-08-16T21:57:37.679460Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# What is going on with those bands?? only appears with voi_lut","metadata":{}},{"cell_type":"markdown","source":"# Is it a matplotlib issue? Try changing the figsize","metadata":{}},{"cell_type":"code","source":"fig, ax = plt.subplots(figsize=(10, 3))\nax.imshow(img3d[:, img3d.shape[1]//2, :])","metadata":{"execution":{"iopub.status.busy":"2021-08-16T21:57:54.764799Z","iopub.execute_input":"2021-08-16T21:57:54.765114Z","iopub.status.idle":"2021-08-16T21:57:54.937304Z","shell.execute_reply.started":"2021-08-16T21:57:54.765087Z","shell.execute_reply":"2021-08-16T21:57:54.936269Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots(figsize=(24, 24))\nax.imshow(img3d[:, img3d.shape[1]//2, :])","metadata":{"execution":{"iopub.status.busy":"2021-08-16T21:58:03.628385Z","iopub.execute_input":"2021-08-16T21:58:03.628884Z","iopub.status.idle":"2021-08-16T21:58:04.124522Z","shell.execute_reply.started":"2021-08-16T21:58:03.628836Z","shell.execute_reply":"2021-08-16T21:58:04.123752Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Does it exist with other colormaps?","metadata":{}},{"cell_type":"code","source":"fig, ax = plt.subplots(figsize=(7, 7))\nax.imshow(img3d[:, img3d.shape[1]//2, :], cmap='gray')","metadata":{"execution":{"iopub.status.busy":"2021-08-16T21:58:25.914347Z","iopub.execute_input":"2021-08-16T21:58:25.914815Z","iopub.status.idle":"2021-08-16T21:58:26.060801Z","shell.execute_reply.started":"2021-08-16T21:58:25.914785Z","shell.execute_reply":"2021-08-16T21:58:26.059998Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}