{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"markdown","source":"

DICOM Processing and Segmentation in Python

\n\n\n\nDICOM is a pain in the neck. It also happens to be very helpful. As clinical radiologists, we expect post-processing, even taking them for granted. However, the magic that occurs behind the scenes is no easy feat, so let’s explore some of that magic.\n\nIn this quest, we will be starting from raw DICOM images. We will extract voxel data from DICOM into numpy arrays, and then perform some low-level operations to normalize and resample the data, made possible using information in the DICOM headers.\n\nThe remainder of the Quest is dedicated to visualizing the data in 1D (by histogram), 2D, and 3D. Finally, we will create segmentation masks that remove all voxel except for the lungs.\n\nProcessing raw DICOM with Python is a little like excavating a dinosaur – you’ll want to have a jackhammer to dig, but also a pickaxe and even a toothbrush for the right situations. Python has all the tools, from pre-packaged imaging process packages handling gigabytes of data at once to byte-level operations on a single voxel.\n\n## Import Packages"},{"metadata":{"trusted":true},"cell_type":"code","source":"# %reload_ext signature\n# %matplotlib inline\n\nimport numpy as np\nimport pydicom\nimport os\nimport matplotlib.pyplot as plt\nfrom glob import glob\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\nimport scipy.ndimage\nfrom skimage import morphology\nfrom skimage import measure\nfrom skimage.transform import resize\nfrom sklearn.cluster import KMeans\nfrom plotly import __version__\nfrom plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot\nimport plotly.figure_factory as ff\nfrom plotly.graph_objs import *\ninit_notebook_mode(connected=True) ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Setting paths\nThen, let's specify a specific DICOM study we can take a closer look. Let's take a look at a chest CT stack from the first patient in the train.csv file given which has some slices that have Pulmonary Embolism (PE). Here we'll use the patient ID 6897fa9de148 from that dataset. \n\n"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"patient_id = '6897fa9de148'\npatient_folder = f'../input/rsna-str-pulmonary-embolism-detection/train/{patient_id}/'\ndata_paths = glob(patient_folder + '/*/*.dcm')\n\n# Print out the first 5 file names to verify we're in the right folder.\nprint (f'Total of {len(data_paths)} DICOM images.\\nFirst 5 filenames:' )\ndata_paths[:5]\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Helper function\n\nHere we make two helper functions.\n\n* `load_scan` will load all DICOM images from a folder into a list for manipulation.\n* The voxel values in the images are raw. `get_pixels_hu` converts raw values into [Houndsfeld units](https://en.wikipedia.org/wiki/Hounsfield_scale)\n\n * The transformation is linear. Therefore, so long as you have a slope and an intercept, you can rescale a voxel value to HU.\n * Both the rescale intercept and rescale slope are stored in the DICOM header at the time of image acquisition (these values are scanner-dependent, so you will need external information)."},{"metadata":{"trusted":true},"cell_type":"code","source":"# \n# Loop over the image files and store everything into a list.\n# \n\ndef load_scan(paths):\n slices = [pydicom.read_file(path ) for path in paths]\n slices.sort(key = lambda x: int(x.InstanceNumber), reverse = True)\n try:\n slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])\n except:\n slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)\n \n for s in slices:\n s.SliceThickness = slice_thickness\n \n return slices\n\ndef get_pixels_hu(scans):\n image = np.stack([s.pixel_array for s in scans])\n # Convert to int16 (from sometimes int16), \n # should be possible as values should always be low enough (<32k)\n image = image.astype(np.int16)\n\n # Set outside-of-scan pixels to 1\n # The intercept is usually -1024, so air is approximately 0\n image[image == -2000] = 0\n \n # Convert to Hounsfield units (HU)\n intercept = scans[0].RescaleIntercept\n slope = scans[0].RescaleSlope\n \n if slope != 1:\n image = slope * image.astype(np.float64)\n image = image.astype(np.int16)\n \n image += np.int16(intercept)\n \n return np.array(image, dtype=np.int16)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Save the loaded images in npy format"},{"metadata":{"trusted":true},"cell_type":"code","source":"patient_id = '0003b3d648eb'\npatient = load_scan(data_paths)\nimgs = get_pixels_hu(patient)\n\n#This is a good time to save the new data set to disk so we don't have to reprocess the stack every time.\nnp.save(f'fullimages_{patient_id}.npy', imgs)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Displaying Images\nThe first thing we should do is to check to see whether the **Houndsfeld Units** are properly scaled and represented. HU's are useful because it is standardized across all CT scans regardless of the absolute number of photons the scanner detector captured. If you need a refresher, here's a quick list of a few useful ones, sourced from Wikipedia.\n\n|Substance|\tHU|\n| :--- | :--- |\n|Air\t|−1000|\n|Lung\t|−500|\n|Fat\t|−100 to −50|\n|Water\t|0|\n|Blood\t|+30 to +70|\n|Muscle\t|+10 to +40|\n|Liver\t|+40 to +60|\n|Bone\t|+700 (cancellous bone) to +3000 (cortical bone)|\n\n\nLet's now create a histogram of all the voxel data in the study."},{"metadata":{"trusted":true},"cell_type":"code","source":"file_used= f'./fullimages_{patient_id}.npy'\nimgs_to_process = np.load(file_used).astype(np.float64) \n\nplt.hist(imgs_to_process.flatten(), bins=100, color='c')\nplt.xlabel(\"Hounsfield Units (HU)\")\nplt.ylabel(\"Frequency\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Critiquing the Histogram\nThe histogram suggests the following:\n\n* There is lots of air (the reason why there are lot of values close to -1000 HU)\n* There is some lung (denoted by the presence of values near -500 HU)\n* There's an abundance of soft tissue, mostly muscle, liver, etc, but there's also some fat. \n* There is only a small bit of bone (seen as a tiny sliver of height between 700-3000)\nThis observation means that we will need to do significant preprocessing if we want to process lesions in the lung tissue because only a tiny bit of the voxels represent lung.\n\nFor reference:\n\n|Substance|\tHU|\n| :--- | :--- |\n|Air\t|−1000|\n|Lung\t|−500|\n|Fat\t|−100 to −50|\n|Water\t|0|\n|Blood\t|+30 to +70|\n|Muscle\t|+10 to +40|\n|Liver\t|+40 to +60|\n|Bone\t|+700 (cancellous bone) to +3000 (cortical bone)|\n\n\nNow if we have some values less than -1000 HU, then we have something strange. Because air really only goes to -1000, so there must be some sort of artifact.\n\nLet's take a look at the actual images.\n\n## Displaying an Image Stack\nWe don't have a lot of screen real estate, so we'll be skipping every 3 slices to get a representative look at the study."},{"metadata":{"trusted":true},"cell_type":"code","source":"file_used= f'./fullimages_{patient_id}.npy'\nimgs_to_process = np.load(file_used).astype(np.float64) \n\ndef sample_stack(stack, rows=6, cols=6, start_with=10, show_every=3):\n fig,ax = plt.subplots(rows,cols,figsize=[18,20])\n for i in range(rows*cols):\n ind = start_with + i*show_every\n ax[int(i/rows),int(i % rows)].set_title(f'slice {ind}')\n ax[int(i/rows),int(i % rows)].imshow(stack[ind],cmap='gray')\n ax[int(i/rows),int(i % rows)].axis('off')\n plt.show()\n\nsample_stack(imgs_to_process)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Resampling\nAlthough we have each individual slices, it is not immediately clear how thick each slice is. Fortunately, this is in the DICOM header."},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f'Slice Thickness: {patient[0].SliceThickness}')\nprint(f'Pixel Spacing (row, col): ({patient[0].PixelSpacing[0]}, {patient[0].PixelSpacing[1]})')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This means we have 2.0 mm slices, and each voxel represents 0.7 mm.\n\nBecause a CT slice is typically reconstructed at 512 x 512 voxels, each slice represents approximately 370 mm of data in length and width.\n\nUsing the metadata from the DICOM we can figure out the size of each voxel as the slice thickness. In order to display the CT in 3D isometric form (which we will do below), and also to compare between different scans, it would be useful to ensure that each slice is resampled in 1x1x1 mm pixels and slices."},{"metadata":{"trusted":true},"cell_type":"code","source":"np.array([float(patient[0].SliceThickness), float(patient[0].PixelSpacing[0]), float(patient[0].PixelSpacing[0])])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"file_used= f'./fullimages_{patient_id}.npy'\nimgs_to_process = np.load(file_used).astype(np.float64)\n\ndef resample(image, scan, new_spacing=[1,1,1]):\n # Determine current pixel spacing\n# spacing = map(float, ([scan[0].SliceThickness] + scan[0].PixelSpacing))\n spacing = np.array([float(patient[0].SliceThickness), \n float(patient[0].PixelSpacing[0]), \n float(patient[0].PixelSpacing[0])])\n\n\n resize_factor = spacing / new_spacing\n new_real_shape = image.shape * resize_factor\n new_shape = np.round(new_real_shape)\n real_resize_factor = new_shape / image.shape\n new_spacing = spacing / real_resize_factor\n \n image = scipy.ndimage.interpolation.zoom(image, real_resize_factor)\n \n return image, new_spacing\n\nprint(f'Shape before resampling: {imgs_to_process.shape}')\nimgs_after_resamp, spacing = resample(imgs_to_process, patient, [1,1,1])\nprint(f'Shape after resampling: {imgs_after_resamp.shape}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def make_mesh(image, threshold=-300, step_size=1):\n\n print(\"Transposing surface\")\n p = image.transpose(2,1,0)\n \n print(\"Calculating surface\")\n verts, faces, norm, val = measure.marching_cubes_lewiner(p, threshold, step_size=step_size, allow_degenerate=True)\n return verts, faces\n\n\ndef plotly_3d(verts, faces):\n x,y,z = zip(*verts) \n \n print(\"Drawing\")\n \n # Make the colormap single color since the axes are positional not intensity. \n# colormap=['rgb(255,105,180)','rgb(255,255,51)','rgb(0,191,255)']\n colormap=['rgb(236, 236, 212)','rgb(236, 236, 212)']\n \n fig = ff.create_trisurf(x=x, y=y, z=z, plot_edges=False,\n colormap=colormap,\n simplices=faces,\n backgroundcolor='rgb(64, 64, 64)',\n title=\"Interactive Visualization\")\n iplot(fig)\n\ndef plt_3d(verts, faces):\n print(\"Drawing\")\n x,y,z = zip(*verts) \n fig = plt.figure(figsize=(10, 10))\n ax = fig.add_subplot(111, projection='3d')\n\n # Fancy indexing: `verts[faces]` to generate a collection of triangles\n mesh = Poly3DCollection(verts[faces], linewidths=0.05, alpha=1)\n face_color = [1, 1, 0.9]\n mesh.set_facecolor(face_color)\n ax.add_collection3d(mesh)\n\n ax.set_xlim(0, max(x))\n ax.set_ylim(0, max(y))\n ax.set_zlim(0, max(z))\n# ax.set_axis_bgcolor((0.7, 0.7, 0.7))\n ax.set_facecolor((0.7,0.7,0.7))\n plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here in the following part, I just have showed a 3D plot of the segmented bone just by using morphological analysis. \n\n* Tune the threshold in the following code block in order to get your desired part segmented. \n* If you want to see an static image of the 3D mesh plot, use the function `plt_3d()`. It takes less time but it is non interactive. That means you can't rotate and twist different axis for better view. \n\n* I you really want to see an interactive 3D mesh plot, just use the function `plotly_3d()` and comment out `plt_3d()`. **Caution: You can only see the interactive 3D plot after you commit the notebook in the kaggle notebook viewer. After the commit even, the interactive 3d plot takes a huge time (~12 mins with my 4Mbps internet)to load. ** "},{"metadata":{"trusted":true},"cell_type":"code","source":"v, f = make_mesh(imgs_after_resamp, threshold = 350) #350 previously default value\nplt_3d(v, f)\n# plotly_3d(v, f)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"## Segmentation\nIf you are interested in chest CTs because you're interested in picking up lung cancers, you're not alone. Machine learning algorithms work a lot better when you can narrowly define what it is looking at. One way to do this is by creating different models for different parts of a chest CT. For instance, a convolutional network for lungs would perform better than a general-purpose network for the whole chest. Therefore, it is often useful to pre-process the image data by auto-detecting the boundaries surrounding a volume of interest.\n\nThe below code will:\n\n* Standardize the pixel value by subtracting the mean and dividing by the standard deviation\n* Identify the proper threshold by creating 2 KMeans clusters comparing centered on soft tissue/bone vs lung/air.\n* Using Erosion) and Dilation) which has the net effect of removing tiny features like pulmonary vessels or noise\n* Identify each distinct region as separate image labels (think the magic wand in Photoshop)\n* Using bounding boxes for each image label to identify which ones represent lung and which ones represent \"every thing else\"\n* Create the masks for lung fields.\n* Apply mask onto the original image to erase voxels outside of the lung fields."},{"metadata":{"trusted":true},"cell_type":"code","source":"def air_removal_mask(dilation):\n labels = measure.label(dilation)\n label_vals = np.unique(labels)\n if labels[0,0] == labels[-1, -1]:\n upper_cut = (labels==labels[0,0])\n mask = np.abs(upper_cut*1 -1) \n else:\n upper_cut = (labels == labels[0,0])\n lower_cut = (labels == labels[-1,-1])\n mask = np.abs((upper_cut + lower_cut )*1 -1) \n return mask\n\n#Standardize the pixel values\ndef make_lungmask(img, display=False):\n row_size= img.shape[0]\n col_size = img.shape[1]\n \n mean = np.mean(img)\n std = np.std(img)\n img = img-mean\n img = img/std\n # Find the average pixel value near the lungs\n # to renormalize washed out images\n middle = img[int(col_size/5):int(col_size/5*4),int(row_size/5):int(row_size/5*4)] \n mean = np.mean(middle) \n max = np.max(img)\n min = np.min(img)\n # To improve threshold finding, I'm moving the \n # underflow and overflow on the pixel spectrum\n img[img==max]=mean\n img[img==min]=mean\n #\n # Using Kmeans to separate foreground (soft tissue / bone) and background (lung/air)\n #\n kmeans = KMeans(n_clusters=2).fit(np.reshape(middle,[np.prod(middle.shape),1]))\n centers = sorted(kmeans.cluster_centers_.flatten())\n threshold = np.mean(centers)\n thresh_img = np.where(imgrow_size/5 and B[2]\n💥 [CT-Scans, DICOM files, Windowing Explained](https://www.kaggle.com/redwankarimsony/ct-scans-dicom-files-windowing-explained)
\n💥 [RSNA-STR-PE [Gradient & Sigmoid Windowing]](https://www.kaggle.com/redwankarimsony/rsna-str-pe-gradient-sigmoid-windowing)
\n💥 [RSNA-STR [✔️3D Stacking ✔️3D Plot ✔️Segmentation]](https://www.kaggle.com/redwankarimsony/rsna-str-3d-stacking-3d-plot-segmentation/edit/run/42517982)
\n💥 [RSNA-STR [DICOM 👉 GIF 👉 npy]](https://www.kaggle.com/redwankarimsony/rsna-str-dicom-gif-npy)
\n💥 [RSNA-STR Pulmonary Embolism [EDA]](https://www.kaggle.com/redwankarimsony/rsna-str-pulmonary-embolism-eda)
\n\n"}],"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":4,"nbformat_minor":4}