{"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>Image Processing - Remove empty spaces, Resize and Save</p></div>\n\n\nDICOM: Digital Imaging and Communications in Medicine\n\nThis notebook reads DICOM images, removes empty spaces in images, resizes images using PIL, and save images to PNG format\n\n","metadata":{}},{"cell_type":"code","source":"!pip install -qU python-gdcm pydicom pylibjpeg","metadata":{"execution":{"iopub.status.busy":"2022-12-02T17:24:13.570518Z","iopub.execute_input":"2022-12-02T17:24:13.571388Z","iopub.status.idle":"2022-12-02T17:24:34.159093Z","shell.execute_reply.started":"2022-12-02T17:24:13.571295Z","shell.execute_reply":"2022-12-02T17:24:34.157818Z"},"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 pandas as pd\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nfrom tqdm.notebook import tqdm\nfrom joblib import Parallel, delayed\nfrom pathlib import Path\n\n#settings\npd.options.display.max_rows = 100\npd.options.display.max_columns = 100\n\n\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n\nimport pytorch_lightning as pl\nrandom_seed=1234\npl.seed_everything(random_seed)","metadata":{"execution":{"iopub.status.busy":"2022-12-02T17:24:34.161546Z","iopub.execute_input":"2022-12-02T17:24:34.161928Z","iopub.status.idle":"2022-12-02T17:24:40.042872Z","shell.execute_reply.started":"2022-12-02T17:24:34.161887Z","shell.execute_reply":"2022-12-02T17:24:40.041754Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df = pd.read_csv(\"/kaggle/input/rsna-breast-cancer-detection/train.csv\")\nprint(train_df.shape)\ntrain_df.head(3)\n","metadata":{"execution":{"iopub.status.busy":"2022-12-02T17:24:40.044645Z","iopub.execute_input":"2022-12-02T17:24:40.04534Z","iopub.status.idle":"2022-12-02T17:24:40.203989Z","shell.execute_reply.started":"2022-12-02T17:24:40.045306Z","shell.execute_reply":"2022-12-02T17:24:40.203151Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def load_img(img_path, resize=True):\n    #img_path=\"/kaggle/input/rsna-breast-cancer-detection/train_images/10038/1967300488.dcm\"\n    dicom = pydicom.dcmread(img_path)\n    img = dicom.pixel_array\n\n    if resize:\n        img = (img - img.min()) / (img.max() - img.min())\n\n        if dicom.PhotometricInterpretation == \"MONOCHROME1\":\n            img = 1 - img\n\n    return img","metadata":{"execution":{"iopub.status.busy":"2022-12-02T17:24:40.206424Z","iopub.execute_input":"2022-12-02T17:24:40.207221Z","iopub.status.idle":"2022-12-02T17:24:40.213273Z","shell.execute_reply.started":"2022-12-02T17:24:40.207187Z","shell.execute_reply":"2022-12-02T17:24:40.212173Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"base_path=\"/kaggle/input/rsna-breast-cancer-detection/train_images\"","metadata":{"execution":{"iopub.status.busy":"2022-12-02T17:24:40.21449Z","iopub.execute_input":"2022-12-02T17:24:40.214839Z","iopub.status.idle":"2022-12-02T17:24:40.226203Z","shell.execute_reply.started":"2022-12-02T17:24:40.214809Z","shell.execute_reply":"2022-12-02T17:24:40.225199Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<p style=\"font-family: 'Mochiy Pop P One';font-size:22px;color:black\" id=1>1. Resize images using PIL package </p>\n<p><a href=#top>back to top</a></p>","metadata":{}},{"cell_type":"code","source":"#Load PIL package\nimport PIL\nfrom PIL import Image\n","metadata":{"execution":{"iopub.status.busy":"2022-12-02T17:24:40.227415Z","iopub.execute_input":"2022-12-02T17:24:40.227957Z","iopub.status.idle":"2022-12-02T17:24:40.24074Z","shell.execute_reply.started":"2022-12-02T17:24:40.227925Z","shell.execute_reply":"2022-12-02T17:24:40.239422Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<p style=\"font-family: 'Mochiy Pop P One';font-size:18px;color:black\" >Define functions to convert DICOM to PIL images  </p>","metadata":{}},{"cell_type":"code","source":"#https://github.com/pydicom/contrib-pydicom/blob/master/viewers/pydicom_PIL.py\ntry:\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":"2022-12-02T17:24:40.244633Z","iopub.execute_input":"2022-12-02T17:24:40.245128Z","iopub.status.idle":"2022-12-02T17:24:40.267455Z","shell.execute_reply.started":"2022-12-02T17:24:40.245083Z","shell.execute_reply":"2022-12-02T17:24:40.266271Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<p style=\"font-family: 'Mochiy Pop P One';font-size:18px;color:black\" >Define functions to remove empty spaces in images  </p>","metadata":{}},{"cell_type":"code","source":"#https://www.kaggle.com/code/jirkaborovec/bloodclots-eda-load-wsi-prune-background?scriptVersionId=101797769\n\ndef prune_image_rows_cols(im, mask, thr=0.990):\n    # delete empty columns\n    for l in reversed(range(im.shape[1])):\n        if (np.sum(mask[:, l]) / float(mask.shape[0])) > thr:\n            im = np.delete(im, l, 1)\n    # delete empty rows\n    for l in reversed(range(im.shape[0])):\n        if (np.sum(mask[l, :]) / float(mask.shape[1])) > thr:\n            im = np.delete(im, l, 0)\n    return im\n\n\ndef mask_median(im, val=255):\n    masks = [None] * 3\n    for c in range(3):\n        masks[c] = im[..., c] >= np.median(im[:, :, c]) - 5\n    mask = np.logical_and(*masks)\n    im[mask, :] = val\n    return im, mask","metadata":{"execution":{"iopub.status.busy":"2022-12-02T17:24:40.268929Z","iopub.execute_input":"2022-12-02T17:24:40.27008Z","iopub.status.idle":"2022-12-02T17:24:40.285854Z","shell.execute_reply.started":"2022-12-02T17:24:40.270034Z","shell.execute_reply":"2022-12-02T17:24:40.284647Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<p style=\"font-family: 'Mochiy Pop P One';font-size:18px;color:black\" >Load a sample image  </p>\n\nThere are two types of `PhotometricInterpretation`: MONOCHROME1; MONOCHROME2\n\n-  MONOCHROME2: the background color is back;\n-  MONOCHROME1: the background color is white;","metadata":{}},{"cell_type":"code","source":"img_path=\"/kaggle/input/rsna-breast-cancer-detection/train_images/10049/1207499426.dcm\"\nimg_path=\"/kaggle/input/rsna-breast-cancer-detection/train_images/10006/1874946579.dcm\"\n#img_path=\"/kaggle/input/rsna-breast-cancer-detection/train_images/10011/270344397.dcm\"\n#img_path=\"/kaggle/input/rsna-breast-cancer-detection/train_images/10130/1672636630.dcm\"\n\ndicom = pydicom.dcmread(img_path)\nimg = dicom.pixel_array\ndicom.PhotometricInterpretation\n#MONOCHROME1; MONOCHROME2","metadata":{"execution":{"iopub.status.busy":"2022-12-02T17:32:15.308314Z","iopub.execute_input":"2022-12-02T17:32:15.308765Z","iopub.status.idle":"2022-12-02T17:32:16.727031Z","shell.execute_reply.started":"2022-12-02T17:32:15.308727Z","shell.execute_reply":"2022-12-02T17:32:16.725876Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"PIL_img=get_PIL_image(dicom)\nimg2 = load_img(img_path, resize=True)\nnp.array(PIL_img).min(), np.array(PIL_img).max(), img2.min(), img2.max()","metadata":{"execution":{"iopub.status.busy":"2022-12-02T17:32:16.729082Z","iopub.execute_input":"2022-12-02T17:32:16.729438Z","iopub.status.idle":"2022-12-02T17:32:18.722837Z","shell.execute_reply.started":"2022-12-02T17:32:16.729406Z","shell.execute_reply":"2022-12-02T17:32:18.721574Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<p style=\"font-family: 'Mochiy Pop P One';font-size:18px;color:black\" >Display and compare different methods of loading images  </p>\n","metadata":{}},{"cell_type":"code","source":"fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(14, 8))\n\n\naxes[0].imshow(img, cmap='gray')\naxes[0].set_title(f'Original Image')\naxes[1].imshow(PIL_img, cmap='gray')\naxes[1].set_title(f'PIL Image')\naxes[2].imshow(img2, cmap='gray')\naxes[2].set_title(f'normalized Image')\n\nplt.suptitle(f\"compare original image and PIL image\", fontsize=15,)\nplt.xticks(fontsize=10)\nplt.yticks(fontsize=10)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-12-02T17:32:18.724593Z","iopub.execute_input":"2022-12-02T17:32:18.724948Z","iopub.status.idle":"2022-12-02T17:32:23.348937Z","shell.execute_reply.started":"2022-12-02T17:32:18.724915Z","shell.execute_reply":"2022-12-02T17:32:23.348112Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<p style=\"font-family: 'Mochiy Pop P One';font-size:18px;color:black\" >Make gray images (number of channels=1) to 3 channels </p>\n\n- Make a numpy array size of width*height to a array of width*height*3\n- For back background images `PhotometricInterpretation==\"MONOCHROME2\"`, convert to white background images","metadata":{}},{"cell_type":"code","source":"if dicom.PhotometricInterpretation==\"MONOCHROME2\": #MONOCHROME1; MONOCHROME2:\n    tmp_img = 255-np.array(PIL_img)\nelse:\n    tmp_img = np.array(PIL_img)\nprint(tmp_img.min(), tmp_img.max())\ntmp_img = np.array([tmp_img, tmp_img, tmp_img]).transpose((1,2,0))\nprint(tmp_img.min(), tmp_img.max())","metadata":{"execution":{"iopub.status.busy":"2022-12-02T17:34:45.496701Z","iopub.execute_input":"2022-12-02T17:34:45.49711Z","iopub.status.idle":"2022-12-02T17:34:45.766364Z","shell.execute_reply.started":"2022-12-02T17:34:45.497076Z","shell.execute_reply":"2022-12-02T17:34:45.765102Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"img6, mask6 = mask_median(tmp_img)\nimg6 = prune_image_rows_cols(img6, mask6)\n\n\nprint(img.size, img6.shape)\nprint(np.asarray(img).min(),np.asarray(img).max())\nprint(img6.min(),img6.max())\n\nfig, axes = plt.subplots(nrows=1, ncols=3, figsize=(16, 8))\naxes[0].imshow(img, cmap='gray')\naxes[0].set_title(f'Original Image')\naxes[1].imshow(PIL_img, cmap='gray')\naxes[1].set_title(f'PIL Image')\naxes[2].imshow(img6, cmap='gray')\naxes[2].set_title(f'.. Image')\n\nplt.suptitle(f\"compare original image and PIL image\", fontsize=15,)\nplt.xticks(fontsize=10)\nplt.yticks(fontsize=10)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-12-02T17:34:47.485792Z","iopub.execute_input":"2022-12-02T17:34:47.486189Z","iopub.status.idle":"2022-12-02T17:36:00.877383Z","shell.execute_reply.started":"2022-12-02T17:34:47.486158Z","shell.execute_reply":"2022-12-02T17:36:00.87606Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<p style=\"font-family: 'Mochiy Pop P One';font-size:18px;color:black\" >Convert image back to channels=1 gray image </p>\n","metadata":{}},{"cell_type":"code","source":"PIL_img = PIL.Image.fromarray(img6[..., 0]) #keep only the first layer","metadata":{"execution":{"iopub.status.busy":"2022-12-02T17:36:00.879819Z","iopub.execute_input":"2022-12-02T17:36:00.880234Z","iopub.status.idle":"2022-12-02T17:36:00.903401Z","shell.execute_reply.started":"2022-12-02T17:36:00.880196Z","shell.execute_reply":"2022-12-02T17:36:00.902009Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#create the thumbnail of the image\n\nif hasattr(Image, 'Resampling'):  # Pillow<8.4.0\n    PIL_img.thumbnail((1024, 1024), resample=Image.Resampling.LANCZOS, reducing_gap=10)\n    if (PIL_img.height< PIL_img.width):\n        PIL_img = PIL_img.transpose(PIL.Image.Transpose.ROTATE_90)\nelse:\n    PIL_img.thumbnail((1024, 1024), resample=Image.LANCZOS, reducing_gap=10)\n    if (PIL_img.height> img.width):\n        PIL_img = PIL_img.transpose(PIL.Image.ROTATE_90)\n        \nnp.array(PIL_img).min(), np.array(PIL_img).max(), img2.min(), img2.max()","metadata":{"execution":{"iopub.status.busy":"2022-12-02T17:36:00.905204Z","iopub.execute_input":"2022-12-02T17:36:00.905732Z","iopub.status.idle":"2022-12-02T17:36:00.992359Z","shell.execute_reply.started":"2022-12-02T17:36:00.905683Z","shell.execute_reply":"2022-12-02T17:36:00.990936Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.array(PIL_img).shape, img.shape, img2.shape","metadata":{"execution":{"iopub.status.busy":"2022-12-02T17:36:00.995572Z","iopub.execute_input":"2022-12-02T17:36:00.996081Z","iopub.status.idle":"2022-12-02T17:36:01.005491Z","shell.execute_reply.started":"2022-12-02T17:36:00.996029Z","shell.execute_reply":"2022-12-02T17:36:01.004091Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f\"PIL image size: {PIL_img.size}, PIL image height: {PIL_img.height}, \\\nPIL image width: {PIL_img.width}, Original image shape: {img.shape}\")\n\nfig, axes = plt.subplots(nrows=1, ncols=3, figsize=(14, 8))\n\n\naxes[0].imshow(img, cmap='gray')\naxes[0].set_title(f'Original Image')\naxes[1].imshow(PIL_img, cmap='gray')\naxes[1].set_title(f'PIL Image')\naxes[2].imshow(img6, cmap='gray')\naxes[2].set_title(f'.. Image')\n\nplt.suptitle(f\"compare original image and PIL image\", fontsize=15,)\nplt.xticks(fontsize=10)\nplt.yticks(fontsize=10)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-12-02T17:36:01.007347Z","iopub.execute_input":"2022-12-02T17:36:01.007879Z","iopub.status.idle":"2022-12-02T17:36:03.207467Z","shell.execute_reply.started":"2022-12-02T17:36:01.007826Z","shell.execute_reply":"2022-12-02T17:36:03.206337Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<p style=\"font-family: 'Mochiy Pop P One';font-size:22px;color:black\" id=2>2. Process all images </p>\n<p><a href=#top>back to top</a></p>","metadata":{}},{"cell_type":"code","source":"target_size = 512\nfor _, row in train_df.iterrows():\n    img_path = f'{base_path}/{row[\"patient_id\"]}/{row[\"image_id\"]}.dcm'\n    if not Path(img_path).exists():\n        print(img_path)\n        continue\n    #img=load_img(img_path, resize=True)\n    dicom = pydicom.dcmread(img_path)\n    PIL_img=get_PIL_image(dicom)\n    \n    \n    #------start to remove empty area-------------------------------------------------------------------\n    if dicom.PhotometricInterpretation==\"MONOCHROME2\": #MONOCHROME1; MONOCHROME2:\n        tmp_img = 255-np.array(PIL_img)\n    else:\n        tmp_img = np.array(PIL_img)\n    tmp_img = np.array([tmp_img, tmp_img, tmp_img]).transpose((1,2,0))\n    img6, mask6 = mask_median(tmp_img)\n    img6 = prune_image_rows_cols(img6, mask6)\n\n    #PIL_img = PIL.Image.fromarray(img6)\n    PIL_img = PIL.Image.fromarray(img6[..., 0]) #keep only the first layer\n\n    #-------end------------------------------------------------------\n    \n    \n    #create the thumbnail of the image\n\n    if hasattr(Image, 'Resampling'):  # Pillow<8.4.0\n        PIL_img.thumbnail((target_size, target_size), resample=Image.Resampling.LANCZOS, reducing_gap=10)\n        if (PIL_img.height< PIL_img.width):\n            PIL_img = PIL_img.transpose(PIL.Image.Transpose.ROTATE_90)\n    else:\n        PIL_img.thumbnail((target_size, target_size), resample=Image.LANCZOS, reducing_gap=10)\n        if (PIL_img.height> img.width):\n            PIL_img = PIL_img.transpose(PIL.Image.ROTATE_90)\n    dest_path = f'{row[\"patient_id\"]}_{row[\"image_id\"]}.png'\n    PIL_img.save(dest_path)","metadata":{"execution":{"iopub.status.busy":"2022-12-02T17:36:04.753512Z","iopub.execute_input":"2022-12-02T17:36:04.754522Z","iopub.status.idle":"2022-12-02T17:37:38.132406Z","shell.execute_reply.started":"2022-12-02T17:36:04.754485Z","shell.execute_reply":"2022-12-02T17:37:38.131283Z"},"trusted":true},"execution_count":null,"outputs":[]}]}