{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Segmenting the right part of the CT image"},{"metadata":{},"cell_type":"markdown","source":"In this short notebook, I want to show you a rather easy way to segement the right part of CT image including just the skull in the middle of the resulting image. This helps your model focus better in the region of interest and also reduce the image size without losing much information. **Please upvote if you found it helpful**"},{"metadata":{},"cell_type":"markdown","source":"I'll be using some code and strategies for correcting bad dicoms files and loading the df from notebooks by [akensert](http://https://www.kaggle.com/akensert/rsna-inceptionv3-keras-tf1-14-0) and [Jeremy Howard](http://https://www.kaggle.com/jhoward/cleaning-the-data-for-rapid-prototyping-fastai). Thank you so much for sharing your great works!"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\nimport cv2\nimport glob\nimport pydicom\nimport numpy as np\nimport pandas as pd\nfrom PIL import Image\nfrom tqdm.notebook import tqdm\nimport matplotlib.pyplot as plt","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"df = pd.read_csv(\"../input/rsna-intracranial-hemorrhage-detection/rsna-intracranial-hemorrhage-detection/stage_2_train.csv\")\ndf[['ID', 'type']] = df['ID'].str.rsplit(\"_\", n=1, expand=True)\ndf.drop_duplicates(['ID', 'type'], inplace=True)\ndf = df.pivot('ID', 'type', 'Label')\ndf.reset_index(inplace=True)\ndf.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Let's take a look at some of the original images**"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_path = \"../input/rsna-intracranial-hemorrhage-detection/rsna-intracranial-hemorrhage-detection/stage_2_train/\"\n\ndef visualize(imgs, axes=None, rows=2, cols=5, **kwargs):\n    if axes is None: _, axes = plt.subplots(rows, cols, **kwargs)\n    for img, ax in zip(imgs, axes.flatten()):\n        ax.imshow(img, cmap='bone')\n        ax.axis(\"off\")\n\nsample_df = df.sample(10, random_state=8)\nsample_df_fnames = [train_path + f\"{id_}.dcm\" for id_ in sample_df['ID'].values]\n\ndcms = [pydicom.dcmread(fname) for fname in sample_df_fnames]\narrs = [dcm.pixel_array for dcm in dcms]\nvisualize(arrs, figsize=(15, 5))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In the next section, I will use the code snippet I found [here](http://https://stackoverflow.com/questions/59865712/segment-photos-of-bio-samples-to-extract-circular-area-of-interest-with-python-i) for a similar problem and I change it a liitle bit to segment the main circular (does not have to be exactly a circle) object in the image. I also return the coordinates and width and height to able to crop the other windows of the CT with the same numbers in order to get fully overlapping images. I will give this function the brain windowed image and use the mentioned numbers to crop other windowed images of the same CT. "},{"metadata":{"trusted":true},"cell_type":"code","source":"def segment_circle(windowed):\n    original = windowed.copy().astype(\"uint8\")\n    mask = np.zeros(original.shape, dtype=np.uint8)\n    gray = original\n    thresh = cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)[1]\n    kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (5,5))\n    close = cv2.morphologyEx(thresh, cv2.MORPH_CLOSE, kernel, iterations=5)\n    cnts = cv2.findContours(close, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)\n    cnts = cnts[0] if len(cnts) == 2 else cnts[1]\n    for c in cnts:\n        peri = cv2.arcLength(c, True)\n        approx = cv2.approxPolyDP(c, 0.04 * peri, True)\n        area = cv2.contourArea(c)\n        if len(approx) > 4 and area > 10000 and area < 500000:\n            ((x, y), r) = cv2.minEnclosingCircle(c)\n            cv2.circle(mask, (int(x), int(y)), int(r), (255, 255, 255), -1)\n            cv2.circle(original, (int(x), int(y)), int(r), (36, 255, 12), 0)\n    x,y,w,h = cv2.boundingRect(mask)\n    mask_ROI = mask[y:y+h, x:x+w]\n    image_ROI = original[y:y+h, x:x+w]\n    result = cv2.bitwise_and(image_ROI, image_ROI, mask=mask_ROI)\n    return result, (x, y, w, h)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I'll use these functions to correct the bad dicoms and window the CT (see the notebooks I put the refernce to at the beginning my notebook)"},{"metadata":{"trusted":true},"cell_type":"code","source":"def correct_dcm(dcm):\n    x = dcm.pixel_array + 1000\n    px_mode = 4096\n    x[x>=px_mode] = x[x>=px_mode] - px_mode\n    dcm.PixelData = x.tobytes()\n    dcm.RescaleIntercept = -1000\n\ndef window_image(img, window_center, window_width):\n    img_min = window_center - window_width // 2\n    img_max = window_center + window_width // 2\n    img = np.clip(img, img_min, img_max)\n    return img","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Here is the main function which reads the dicoms, windows them and then segments accordingly**"},{"metadata":{"trusted":true},"cell_type":"code","source":"def preprocess_and_segment(path, norm='min_max', make_correct=True):    \n    dcm = pydicom.dcmread(path)\n    if make_correct and (dcm.BitsStored == 12) and (dcm.PixelRepresentation == 0) and (int(dcm.RescaleIntercept) > -100):\n        correct_dcm(dcm)\n\n    sample = dcm.pixel_array.astype(\"float32\") * dcm.RescaleSlope + dcm.RescaleIntercept\n    \n    brain_window = window_image(sample, 40, 80)\n    dural_window = window_image(sample, 80, 200)\n    bone_window = window_image(sample, 600, 2800)\n    \n    brain_window, (x,y,w,h) = segment_circle(brain_window)\n    if brain_window is None:\n        return np.zeros((512, 512, 3))\n    dural_window = dural_window[y:y+h, x:x+w]\n    bone_window = bone_window[y:y+h, x:x+w]\n    \n    # Just some normalization. You can ignore this\n    if norm == \"min_max\":\n        brain_window = (brain_window - (0.)) / 80.\n        dural_window = (dural_window - (-20.)) / 200.\n        bone_window = (bone_window - (-1200.)) / 2800.\n        img_3ch = np.dstack([brain_window, dural_window, bone_window]).astype(\"float32\")\n    elif norm == \"hu\":\n        img_3ch = np.dstack([brain_window, dural_window, bone_window]).astype(\"float32\") / 4096.\n    elif norm == 'none':\n        img_3ch = np.dstack([brain_window, dural_window, bone_window]).astype(\"float32\")\n    \n    return img_3ch","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**The original images (again, for comparison)**"},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_df = df.sample(10, random_state=8)\nsample_df_fnames = [train_path + f\"{id_}.dcm\" for id_ in sample_df['ID'].values]\n\ndcms = [pydicom.dcmread(fname) for fname in sample_df_fnames]\narrs = [dcm.pixel_array for dcm in dcms]\nvisualize(arrs, figsize=(15, 5))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**After segmenting (cropping) the right part of the image**"},{"metadata":{"trusted":true},"cell_type":"code","source":"proccessed = [preprocess_and_segment(fname, norm=\"min_max\") for fname in sample_df_fnames]\nvisualize([img[..., 0] for img in proccessed], figsize=(15, 5))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As you can see, the segmentation function successfully cropped the right part of the image and the skull is in the middle of the resulting images. When the brain tissue is really low or in some other rare situations, the segmentation function could fail (which is not frequent and happens mostly when the image label is 0 with no hemorrhage) and returns **None** which in this case I'm retuning a black image (np array full of zeros). You can see an example where the function failed to segment in the first image of second row."},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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}