{"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":"

Understanding DICOM (And how to use pydicom)

\n\n

DICOM stands for Digital Imaging and Communications in Medicine. Pretty straightforward that all data generated by machines(X-rays, CT-scans, Radiotherapy images etc.) in healthcare are stored in this format.

\n\n
\n\n

Usually a DICOM file or .dcm file will contain:

\n\n*

File meta-data (Yes, as it's content! Pretty surprising for cybersec brethren innit)

\n*

Information about the patient

\n*

Times of procedures/scans, patient-positions etc.

\n*

Photometric interpretation (more on this later)

\n*

Pixel data (or simply the medical image which is usually a set of monochrome channels)

\n\n
","metadata":{}},{"cell_type":"code","source":"!pip install -Uqq pydicom\nimport pydicom\n\nimport cv2\nimport os\nfrom tqdm.notebook import tqdm\nfrom pathlib import Path\n\nfrom PIL import Image\nfrom matplotlib import pyplot as plt\nimport numpy as np","metadata":{"execution":{"iopub.status.busy":"2023-08-21T06:46:44.894523Z","iopub.execute_input":"2023-08-21T06:46:44.895018Z","iopub.status.idle":"2023-08-21T06:47:02.786145Z","shell.execute_reply.started":"2023-08-21T06:46:44.894979Z","shell.execute_reply":"2023-08-21T06:47:02.784725Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n\n

Loading, Accessing & Understanding a .dcm file OR pydicom DataElement

\n\n
\n\n

Reference(s): \n\n

1. Pydicom Documentation: DataElement

\n

2. pydicom docs: dataset basics

","metadata":{}},{"cell_type":"code","source":"path = Path('/kaggle/input/rsna-2023-abdominal-trauma-detection/train_images/10005/18667/100.dcm') # using a random file from the training set\ndcm = pydicom.dcmread(path) # read the file\nlen(dcm) ","metadata":{"execution":{"iopub.status.busy":"2023-08-21T06:47:02.788584Z","iopub.execute_input":"2023-08-21T06:47:02.789064Z","iopub.status.idle":"2023-08-21T06:47:02.835473Z","shell.execute_reply.started":"2023-08-21T06:47:02.789017Z","shell.execute_reply":"2023-08-21T06:47:02.834232Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**🦄 Ooof... has 29 objects inside! Let's see each one**","metadata":{}},{"cell_type":"code","source":"print(\"S.No | Name | VR | VM | Value\\n-----------------------------\")\nfor i, itm in enumerate(dcm):\n if itm.name == \"Pixel Data\":\n print(i+1, itm.name, itm.VR, itm.VM, \"| Pixel Data Not printed due to huge size\")\n else:\n print(i+1, itm.name, itm.VR, itm.VM, itm.value)","metadata":{"execution":{"iopub.status.busy":"2023-08-21T06:47:02.837045Z","iopub.execute_input":"2023-08-21T06:47:02.837414Z","iopub.status.idle":"2023-08-21T06:47:02.846455Z","shell.execute_reply.started":"2023-08-21T06:47:02.837386Z","shell.execute_reply":"2023-08-21T06:47:02.845464Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n\n

Members of each item present in the pydicom dataElement

\n\n

1. tag: the element’s tag (A unique integer for different things like patient_name) (Can be used as a key to access the corresponding item)

\n

2. name: the element's name (or the parsed tag) (Can also be used as a key)

\n

3. VR: Value Representations → a two letter str that describes the format of the stored value.\n

4. VM: the element’s Value Multiplicity as an int\n

5. value: the element’s actual value \n\n

\n\n

Medical Images == Even More Fun Stuff

\n\n

Can directly call pixel_array on pydicom DataElement!

\n\n
","metadata":{}},{"cell_type":"code","source":"medImg = dcm.pixel_array # Can treat dcm like a class object with member variables\nmedImg.shape, medImg.dtype, type(medImg), np.amax(medImg), np.amin(medImg)","metadata":{"execution":{"iopub.status.busy":"2023-08-21T06:47:02.849373Z","iopub.execute_input":"2023-08-21T06:47:02.850564Z","iopub.status.idle":"2023-08-21T06:47:02.878558Z","shell.execute_reply.started":"2023-08-21T06:47:02.850522Z","shell.execute_reply":"2023-08-21T06:47:02.877498Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**🦄 Whoa 16 bit image with negative min value!! Gotta fix this and store PNGs for our models**","metadata":{}},{"cell_type":"code","source":"data = medImg\n# Normalize to 0, 1 range\ndata = (data - data.min()) / (data.max() - data.min())","metadata":{"execution":{"iopub.status.busy":"2023-08-21T06:47:02.880495Z","iopub.execute_input":"2023-08-21T06:47:02.880833Z","iopub.status.idle":"2023-08-21T06:47:02.891617Z","shell.execute_reply.started":"2023-08-21T06:47:02.880806Z","shell.execute_reply":"2023-08-21T06:47:02.890617Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n\n`MONOCHROME1` indicates that the grayscale ranges from bright to dark with ascending pixel values \n\nwhereas, \n\n`MONOCHROME2` ranges from dark to bright with ascending pixel values.\n\n
","metadata":{}},{"cell_type":"code","source":"if dcm.PhotometricInterpretation == \"MONOCHROME1\":\n data = 1 - data # invert to view in matplotlib\nelif dcm.PhotometricInterpretation == \"MONOCHROME2\":\n data = data\n\ndata = (data * 255).astype(np.uint8) # Make it a standard 8-bit array | for PNG/JPG format\n\nplt.imshow(data, cmap=plt.cm.bone)\nplt.axis('off')\nplt.colorbar()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-08-21T06:47:02.893293Z","iopub.execute_input":"2023-08-21T06:47:02.893612Z","iopub.status.idle":"2023-08-21T06:47:03.409424Z","shell.execute_reply.started":"2023-08-21T06:47:02.893586Z","shell.execute_reply":"2023-08-21T06:47:03.408334Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n\n

If you made it till here, here's a general utility function for manipulating & saving DICOM Images for your Machine Learning Models!

\n\n
","metadata":{}},{"cell_type":"code","source":"def preprocess_dicom(dcm, isSave=False, savePath=None):\n # NOTE: Assumes opencv is imported as cv2 & numpy as np\n if isSave:\n assert savePath is not None, \"Please provide a save path if you want to save the preprocessed image\"\n \n data = dcm.pixel_array\n\n # normalize to fix negative values & also the standard pixel input to any model is [0,1]\n data = (data - data.min()) / (data.max() - data.min())\n\n if dcm.PhotometricInterpretation == \"MONOCHROME1\":\n data = 1 - data # invert to correct\n elif dcm.PhotometricInterpretation == \"MONOCHROME2\":\n data = data\n\n if isSave:\n data = (data * 255).astype(np.uint8)\n cv2.imwrite(savePath + str(dcm.PatientID) + \".png\", data)\n\n return data","metadata":{"execution":{"iopub.status.busy":"2023-08-21T06:47:03.410587Z","iopub.execute_input":"2023-08-21T06:47:03.41094Z","iopub.status.idle":"2023-08-21T06:47:03.425182Z","shell.execute_reply.started":"2023-08-21T06:47:03.410909Z","shell.execute_reply":"2023-08-21T06:47:03.424157Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.imshow(preprocess_dicom(dcm))\nplt.colorbar()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-08-21T06:47:03.426143Z","iopub.execute_input":"2023-08-21T06:47:03.426457Z","iopub.status.idle":"2023-08-21T06:47:03.867625Z","shell.execute_reply.started":"2023-08-21T06:47:03.426431Z","shell.execute_reply":"2023-08-21T06:47:03.866393Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Plotting The 3D CT-scans \n\n**[Check this awesome answer about how to get meaningful 3D volumes from DICOM images](https://stackoverflow.com/questions/8756096/window-width-and-center-calculation-of-dicom-image/8765366#8765366)**\n\n> Load multiple images & build 3d image $\\rightarrow$ Plot 3D Volume $\\rightarrow$ Re-slice for angular views","metadata":{}},{"cell_type":"code","source":"from tqdm.notebook import tqdm","metadata":{"execution":{"iopub.status.busy":"2023-08-21T06:47:03.869084Z","iopub.execute_input":"2023-08-21T06:47:03.869515Z","iopub.status.idle":"2023-08-21T06:47:03.873862Z","shell.execute_reply.started":"2023-08-21T06:47:03.869485Z","shell.execute_reply":"2023-08-21T06:47:03.872812Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Load multiple images & Build 3D Volume","metadata":{}},{"cell_type":"code","source":"SAMPLE_PATH = \"/kaggle/input/rsna-2023-abdominal-trauma-detection/train_images/10005/18667\"","metadata":{"execution":{"iopub.status.busy":"2023-08-21T06:47:03.877002Z","iopub.execute_input":"2023-08-21T06:47:03.87734Z","iopub.status.idle":"2023-08-21T06:47:03.888125Z","shell.execute_reply.started":"2023-08-21T06:47:03.877311Z","shell.execute_reply":"2023-08-21T06:47:03.887101Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dcm_fnames = []\nskip_fnames = []\nvol_3d = []\n\nfor fname in tqdm(os.listdir(SAMPLE_PATH)):\n dcm = pydicom.dcmread(os.path.join(SAMPLE_PATH,fname))\n if hasattr(dcm, \"RescaleIntercept\") and hasattr(dcm, \"RescaleSlope\"):\n # Rescale params\n intercept = float(dcm.RescaleIntercept)\n slope = float(dcm.RescaleSlope)\n \n # Clipping params\n center = int(dcm.WindowCenter)\n width = int(dcm.WindowWidth)\n low = center - width / 2\n high = center + width / 2 \n \n image = np.clip(dcm.pixel_array * slope + intercept, low, high)\n \n vol_3d.append(image)\n \n dcm_fnames.append(fname)\n else:\n skip_fnames.append(fname)\n\nvolume = np.stack(vol_3d, axis=0)\n\nprint(f\"Skipped {len(skip_fnames)} samples\\nUsed: {len(dcm_fnames)}\\n\")\nprint(\"Contructed Volume.shape:\", volume.shape)","metadata":{"execution":{"iopub.status.busy":"2023-08-21T06:47:03.889641Z","iopub.execute_input":"2023-08-21T06:47:03.889989Z","iopub.status.idle":"2023-08-21T06:47:07.111281Z","shell.execute_reply.started":"2023-08-21T06:47:03.889935Z","shell.execute_reply":"2023-08-21T06:47:07.109852Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Show The 3D Volume\n\n\n#### [StackOverflow Reference for 3D-Visualization](https://stackoverflow.com/questions/56035562/3d-dicom-visualisation-in-python) (OUTDATED) $\\rightarrow$ Updated using [skimage doc](https://scikit-image.org/docs/dev/api/skimage.measure.html#skimage.measure.marching_cubes) for marching cubes","metadata":{}},{"cell_type":"code","source":"from mpl_toolkits.mplot3d.art3d import Poly3DCollection\nimport numpy as np\nfrom skimage import measure\n\n# WARNING: This method can take a long time!\ndef plot_3d(image): \n p = image#.transpose(2,1,0)\n verts, faces, normals, values = measure.marching_cubes(p, \n step_size=2, \n# level=-50, \n method='lewiner',\n allow_degenerate=True)\n fig = plt.figure(figsize=(10, 10))\n ax = fig.add_subplot(111, projection='3d')\n mesh = Poly3DCollection(verts[faces], alpha=0.1)\n face_color = [0.5, 0.5, 1]\n mesh.set_facecolor(face_color)\n ax.add_collection3d(mesh)\n ax.set_xlim(0, p.shape[0])\n ax.set_ylim(0, p.shape[1])\n ax.set_zlim(0, p.shape[2])\n\n plt.show()\n\n# plot_3d(volume)","metadata":{"execution":{"iopub.status.busy":"2023-08-21T06:47:07.113056Z","iopub.execute_input":"2023-08-21T06:47:07.113585Z","iopub.status.idle":"2023-08-21T06:47:07.671505Z","shell.execute_reply.started":"2023-08-21T06:47:07.113543Z","shell.execute_reply":"2023-08-21T06:47:07.670461Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Reslice & Display with different angles (Sagittal, Coronal & Axial) $\\rightarrow$ Reference Code: [Awesome NB by Parham Mostame](https://www.kaggle.com/code/parhammostame/construct-3d-arrays-from-dcm-nii-3-view-angles) \n","metadata":{}},{"cell_type":"code","source":"def plot_image_with_seg(volume, volume_seg=[], orientation='Coronal', num_subplots=20):\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) \n \n \nplot_image_with_seg(volume, orientation='Coronal', num_subplots=10)\nplot_image_with_seg(volume, orientation='Sagittal', num_subplots=10)\nplot_image_with_seg(volume, orientation='Axial', num_subplots=10)","metadata":{"execution":{"iopub.status.busy":"2023-08-21T06:47:26.961935Z","iopub.execute_input":"2023-08-21T06:47:26.962378Z","iopub.status.idle":"2023-08-21T06:47:29.641411Z","shell.execute_reply.started":"2023-08-21T06:47:26.962345Z","shell.execute_reply":"2023-08-21T06:47:29.640171Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Thanks for checking this NB out ❤️ ","metadata":{}}]}