{"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":"<div align=\"center\"><p style=\"font-family: 'Mochiy Pop P One';font-size:32px;color:black\" id=top>DICOM image processing</p></div>\n\n\nDICOM: Digital Imaging and Communications in Medicine\n\nThis notebook reads DICOM files and crops them to exclude the background. \nIt also removes orientation labels that are common in DICOM images, so that the labels do not affect cropping the image.\nDeveloped using the RSNA Screening Mammography Breast Cancer Detection dataset.\n\nForked from XXXXYYYY80008.\n\n","metadata":{}},{"cell_type":"code","source":"!pip install -qU python-gdcm pydicom pylibjpeg","metadata":{"execution":{"iopub.status.busy":"2023-01-26T22:45:51.300376Z","iopub.execute_input":"2023-01-26T22:45:51.30113Z","iopub.status.idle":"2023-01-26T22:46:00.722957Z","shell.execute_reply.started":"2023-01-26T22:45:51.301024Z","shell.execute_reply":"2023-01-26T22:46:00.72104Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#!pip install Pillow","metadata":{"execution":{"iopub.status.busy":"2023-01-26T22:46:00.725342Z","iopub.execute_input":"2023-01-26T22:46:00.725684Z","iopub.status.idle":"2023-01-26T22:46:00.731126Z","shell.execute_reply.started":"2023-01-26T22:46:00.725656Z","shell.execute_reply":"2023-01-26T22:46:00.729872Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport cv2\nimport glob\nimport gdcm\nimport copy\nimport pydicom\nimport numpy as np\nimport matplotlib.pyplot as plt\n\nfrom pathlib import Path\n\nimport PIL\nfrom PIL import Image\nfrom PIL import ImageOps\nfrom PIL import ImageFilter\n\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n","metadata":{"execution":{"iopub.status.busy":"2023-01-26T22:46:00.732725Z","iopub.execute_input":"2023-01-26T22:46:00.733225Z","iopub.status.idle":"2023-01-26T22:46:00.894714Z","shell.execute_reply.started":"2023-01-26T22:46:00.733179Z","shell.execute_reply":"2023-01-26T22:46:00.892907Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The following cell in from pydicom_PIL.py.  \nhttps://github.com/pydicom/contrib-pydicom/blob/master/viewers/pydicom_PIL.py\n","metadata":{}},{"cell_type":"code","source":"try:\n    import PIL.Image\n    have_PIL = True\nexcept ImportError:\n    have_PIL = False\n\ntry:\n    import numpy as np\n    have_numpy = True\nexcept ImportError:\n    have_numpy = False\n\n\ndef get_LUT_value(data, window, level):\n    \"\"\"Apply the RGB Look-Up Table for the given\n       data and window/level value.\"\"\"\n    if not have_numpy:\n        raise ImportError(\"Numpy is not available.\"\n                          \"See http://numpy.scipy.org/\"\n                          \"to download and install\")\n\n    return np.piecewise(data,\n                        [data <= (level - 0.5 - (window - 1) / 2),\n                         data > (level - 0.5 + (window - 1) / 2)],\n                        [0, 255, lambda data: ((data - (level - 0.5)) /\n                         (window - 1) + 0.5) * (255 - 0)])\n\n\ndef get_PIL_image(dataset):\n    \"\"\"Get Image object from Python Imaging Library(PIL)\"\"\"\n    if not have_PIL:\n        raise ImportError(\"Python Imaging Library is not available. \"\n                          \"See http://www.pythonware.com/products/pil/ \"\n                          \"to download and install\")\n\n    if ('PixelData' not in dataset):\n        raise TypeError(\"Cannot show image -- DICOM dataset does not have \"\n                        \"pixel data\")\n    # can only apply LUT if these window info exists\n    if ('WindowWidth' not in dataset) or ('WindowCenter' not in dataset):\n        bits = dataset.BitsAllocated\n        samples = dataset.SamplesPerPixel\n        if bits == 8 and samples == 1:\n            mode = \"L\"\n        elif bits == 8 and samples == 3:\n            mode = \"RGB\"\n        elif bits == 16:\n            # not sure about this -- PIL source says is 'experimental'\n            # and no documentation. Also, should bytes swap depending\n            # on endian of file and system??\n            mode = \"I;16\"\n        else:\n            raise TypeError(\"Don't know PIL mode for %d BitsAllocated \"\n                            \"and %d SamplesPerPixel\" % (bits, samples))\n\n        # PIL size = (width, height)\n        size = (dataset.Columns, dataset.Rows)\n\n        # Recommended to specify all details\n        # by http://www.pythonware.com/library/pil/handbook/image.htm\n        im = PIL.Image.frombuffer(mode, size, dataset.PixelData,\n                                  \"raw\", mode, 0, 1)\n\n    else:\n        ew = dataset['WindowWidth']\n        ec = dataset['WindowCenter']\n        ww = int(ew.value[0] if ew.VM > 1 else ew.value)\n        wc = int(ec.value[0] if ec.VM > 1 else ec.value)\n        image = get_LUT_value(dataset.pixel_array, ww, wc)\n        # Convert mode to L since LUT has only 256 values:\n        #   http://www.pythonware.com/library/pil/handbook/image.htm\n        im = PIL.Image.fromarray(image).convert('L')\n\n    return im\n\n\ndef show_PIL(dataset):\n    \"\"\"Display an image using the Python Imaging Library (PIL)\"\"\"\n    im = get_PIL_image(dataset)\n    im.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-26T22:46:00.898071Z","iopub.execute_input":"2023-01-26T22:46:00.89846Z","iopub.status.idle":"2023-01-26T22:46:00.913307Z","shell.execute_reply.started":"2023-01-26T22:46:00.898431Z","shell.execute_reply":"2023-01-26T22:46:00.912282Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The following routines crop images, and by default remove labels.  \nThe image cropping function is 16x faster than the one in the notebook that I forked that used array manipulation.  ","metadata":{}},{"cell_type":"code","source":"def crop_median_image(img, val=255):\n    \"\"\"Masks and crops an image based on the median value.\n    \n    PARAMETERS:\n        img - An Image\n    RETURNS:\n        A cropped Image.\n    \"\"\"\n    \n    img_ary = np.array(img)\n    masks = [None] * 3\n    for c in range(3):\n        masks[c] = img_ary[..., c] >= np.median(img_ary[:, :, c]) - 5\n    mask = np.logical_and(*masks)\n    img_ary[mask, :] = val\n    masked_img = Image.fromarray(img_ary)\n\n    mask = np.array(mask, dtype=np.uint8)\n    # smooth the mask\n    mask_img = Image.fromarray(mask)\n    mask_img = mask_img.filter(ImageFilter.GaussianBlur(radius=255))\n    \n    mask = np.array(np.array(mask_img) > 0, dtype=np.bool8)\n    \n    mask_inv = ImageOps.invert(Image.fromarray(mask))\n    mask_bb = mask_inv.getbbox()\n    cropped = masked_img.crop(mask_bb)\n    return cropped, mask_inv\n\ndef dicom_to_gray_image(dicom):\n    \"\"\"Reads a DICOM image and returns a grayscale Image.\n    \n    PARAMETERS:\n        dicom - A DICOM image.\n    RETURNS:\n        An Image.\n    \"\"\"\n    \n    img=get_PIL_image(dicom)\n    \n    if dicom.PhotometricInterpretation==\"MONOCHROME2\": #MONOCHROME1; MONOCHROME2:\n        tmp_img = 255-np.array(img)\n    else:\n        tmp_img = np.array(img)\n    tmp_img = np.array([tmp_img, tmp_img, tmp_img]).transpose((1,2,0))\n    return Image.fromarray(tmp_img)\n                           \ndef dicom_to_cropped_image(dicom, want_mask=False, keep_labels=False):\n    \"\"\"Convert a dicom to a cropped Image with a white background.\n    \n    PARAMETER:\n        dicom - A dicom image.\n        want_mask - If true, returns the cropped image and the mask image as a tuple. Default is False.\n        keep_labels - If true, does not remove labels from images. Default is False.\n    RETURNS:\n        A cropped Image, or, if want_mask is True, the cropped Image and Mask as a tuple.\n    \"\"\"\n\n    tmp_img = dicom_to_gray_image(dicom)\n\n    cropped_img, mask = crop_median_image(tmp_img)\n    cropped_ary = np.array(cropped_img)\n   \n    cropped_img = PIL.Image.fromarray(cropped_ary[..., 0]) #keep only the first layer\n\n    #create the thumbnail of the image\n\n    if hasattr(Image, 'Resampling'):  # Pillow<8.4.0\n        cropped_img.thumbnail((1024, 1024), resample=Image.Resampling.LANCZOS, reducing_gap=10)\n        if (cropped_img.height< cropped_img.width):\n            cropped_img = cropped_img.transpose(PIL.Image.Transpose.ROTATE_90)\n    else:\n        cropped_img.thumbnail((1024, 1024), resample=Image.LANCZOS, reducing_gap=10)\n        if (cropped_img.height> img.width):\n            cropped_img = cropped_img.transpose(PIL.Image.ROTATE_90)\n            \n    if want_mask:\n        return(cropped_img, mask)\n\n    return cropped_img","metadata":{"execution":{"iopub.status.busy":"2023-01-26T22:46:00.915595Z","iopub.execute_input":"2023-01-26T22:46:00.916365Z","iopub.status.idle":"2023-01-26T22:46:01.1419Z","shell.execute_reply.started":"2023-01-26T22:46:00.916317Z","shell.execute_reply":"2023-01-26T22:46:01.141087Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Test Code\nCode to test the above.","metadata":{"execution":{"iopub.status.busy":"2023-01-26T17:51:36.538405Z","iopub.execute_input":"2023-01-26T17:51:36.538829Z","iopub.status.idle":"2023-01-26T17:51:36.545565Z","shell.execute_reply.started":"2023-01-26T17:51:36.53879Z","shell.execute_reply":"2023-01-26T17:51:36.544225Z"}}},{"cell_type":"code","source":"#img_paths = ['/kaggle/input/rsna-breast-cancer-detection/train_images/34867/1816992952.dcm']\nimg_paths = glob.glob(\"/kaggle/input/rsna-breast-cancer-detection/train_images/**/*.dcm\")[:10]\n\nfor image_path in img_paths:\n    print(image_path)\n    dicom = pydicom.dcmread(image_path)\n    original_image = dicom_to_gray_image(dicom)\n    image, mask = dicom_to_cropped_image(dicom, want_mask=True)\n\n    print(original_image.size, image.size)\n\n    fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(16, 8))\n    axes[0].imshow(original_image, cmap='gray')\n    axes[0].set_title(f'Original image {original_image.size}')\n    axes[1].imshow(mask, cmap='gray')\n    axes[1].set_title(f'the mask')\n    axes[2].imshow(image, cmap='gray')\n    axes[2].set_title(f'Cropped image {image.size}')\n    plt.show()\n    \n    # save the image, if desired\n#     base_name = os.path.splitext(image_path)[0]\n#     image.save(base_name + \".png\")\n","metadata":{"execution":{"iopub.status.busy":"2023-01-26T22:46:01.143097Z","iopub.execute_input":"2023-01-26T22:46:01.143536Z","iopub.status.idle":"2023-01-26T22:46:57.407167Z","shell.execute_reply.started":"2023-01-26T22:46:01.143507Z","shell.execute_reply":"2023-01-26T22:46:57.406051Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"image_path","metadata":{"execution":{"iopub.status.busy":"2023-01-26T22:46:57.408705Z","iopub.execute_input":"2023-01-26T22:46:57.409152Z","iopub.status.idle":"2023-01-26T22:46:57.416367Z","shell.execute_reply.started":"2023-01-26T22:46:57.409116Z","shell.execute_reply":"2023-01-26T22:46:57.415051Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}