{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Introduction "},{"metadata":{},"cell_type":"markdown","source":"In this kernel notebook I show how to setup a DicomImage Databunch object with visualization which can be used with the fastai v1 library."},{"metadata":{},"cell_type":"markdown","source":"# Setup"},{"metadata":{},"cell_type":"markdown","source":"We start with the import of the libraries we need."},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:06.768779Z","start_time":"2019-10-18T13:41:05.627063Z"},"trusted":true},"cell_type":"code","source":"import pandas as pd\nimport pydicom\nfrom fastai.vision import *\nfrom fastai.data_block import _maybe_squeeze\nfrom collections import defaultdict","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"You can use the [turbo colormap mentioned in a kernel notebook discussion](https://www.kaggle.com/jhoward/don-t-see-like-a-radiologist-fastai) by downloading the turbo colormap from the [GitHub gist](https://gist.github.com/FedeMiorelli/640bbc66b2038a14802729e609abfe89) (or via your shell with `wget https://gist.githubusercontent.com/FedeMiorelli/640bbc66b2038a14802729e609abfe89/raw/34a26667e1528c9e4465cbc0be30d10cbe8d4a40/turbo_colormap_mpl.py`). *Unfortunately, I cannot turn on the Internet in this kernel, so this is not working.*"},{"metadata":{"trusted":true},"cell_type":"code","source":"#!wget \"https://gist.githubusercontent.com/FedeMiorelli/640bbc66b2038a14802729e609abfe89/raw/34a26667e1528c9e4465cbc0be30d10cbe8d4a40/turbo_colormap_mpl.py\"","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:06.77368Z","start_time":"2019-10-18T13:41:06.769827Z"},"trusted":true},"cell_type":"code","source":"#import turbo_colormap_mpl","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:06.778038Z","start_time":"2019-10-18T13:41:06.774628Z"},"trusted":true},"cell_type":"code","source":"# show fastai version\n__version__","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Data"},{"metadata":{},"cell_type":"markdown","source":"## Preparation"},{"metadata":{},"cell_type":"markdown","source":"The data preparation is heavily based on [Radeks fastai starter kit](https://www.kaggle.com/c/rsna-intracranial-hemorrhage-detection/discussion/109649#latest-651858) from his [GitHub repo](https://github.com/radekosmulski/rsna-intracranial). (Thanks Radek!)"},{"metadata":{"trusted":true},"cell_type":"code","source":"df = pd.read_csv('../input/rsna-intracranial-hemorrhage-detection/stage_1_train.csv')\ndf['fn'] = df.ID.apply(lambda x: '_'.join(x.split('_')[:2]) + '.dcm')\ndf.columns = ['ID', 'probability', 'fn']\ndf['label'] = df.ID.apply(lambda x: x.split('_')[-1])\ndf.drop_duplicates('ID', inplace=True)\npivot = df.pivot(index='fn', columns='label', values='probability')\npivot.reset_index(inplace=True)\nd = defaultdict(list)\n\nfor fn in df.fn.unique(): d[fn]\n\nfor tup in df.itertuples():\n    if tup.probability: d[tup.fn].append(tup.label)\n        \nks, vs = [], []\n\nfor k, v in d.items():\n    ks.append(k), vs.append(' '.join(v))\n    \ndf_train = pd.DataFrame(data={'fn': ks, 'labels': vs})","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:07.03977Z","start_time":"2019-10-18T13:41:07.014681Z"},"trusted":true},"cell_type":"code","source":"df_train.labels.fillna('', inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:07.047363Z","start_time":"2019-10-18T13:41:07.040657Z"},"trusted":true},"cell_type":"code","source":"# safety check of df shape\nassert df_train.shape == (674258, 2)","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:07.085556Z","start_time":"2019-10-18T13:41:07.048347Z"},"trusted":true},"cell_type":"code","source":"# remove damaged dicom file from df\ndf_train = df_train[df_train['fn'] != 'ID_6431af929.dcm']","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:07.088761Z","start_time":"2019-10-18T13:41:07.087017Z"},"trusted":true},"cell_type":"code","source":"# safety check of df shape\nassert df_train.shape == (674258-1, 2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Dicom classes & functions"},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:07.093607Z","start_time":"2019-10-18T13:41:07.089799Z"},"trusted":true},"cell_type":"code","source":"# change here the colormap for the entire notebook\n#cmap='turbo'\ncmap='jet'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We setup a a DicomImage class to handle our dicom images, including a open function and a show function. We setup a `open_dicom_image` function to generate images with 3 channels from the image data with 1 channel by expanding. With that setup we can use pretrained networks based on 3 channel/RGB images (this can be disabled in `open_dicom_image` by setting `expand=False`.) In the setup outlined we clamp the rescaled dicom image data values to be between -1024 and 1024 (this can be changed in `open_dicom_image` by changing `clamp_min` and `clamp_max`)."},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:07.104654Z","start_time":"2019-10-18T13:41:07.099872Z"},"trusted":true},"cell_type":"code","source":"class DicomImage(Image):\n    \"DicomImage to support applying transforms to image data in `px`.\"\n    def __init__(self, px:Tensor):\n        self._px = px\n        self._logit_px=None\n        self._flow=None\n        self._affine_mat=None\n        self.sample_kwargs = {}\n    \n    def _repr_image_format(self, format_str):\n        with BytesIO() as str_buffer:\n            plt.imsave(str_buffer, image2np(self.px[0].unsqueeze(0)), # We show only one channel!\n                       format=format_str, cmap=cmap)\n            return str_buffer.getvalue()\n        \n    def clone(self):\n        \"Mimic the behavior of torch.clone for `Image` objects.\"\n        return self.__class__(self.px.clone())\n\n    @property\n    def data(self)->TensorImage:\n        \"Return this images pixels as a tensor.\"\n        return self.px\n    \n    def show(self, ax:plt.Axes=None, figsize:tuple=(3,3), title:Optional[str]=None, hide_axis:bool=True,\n              cmap:str=cmap, y:Any=None, **kwargs):\n        \"Show image on `ax` with `title`, overlaid with optional `y`\"\n        ax = show_dicom_image(self, ax=ax, hide_axis=hide_axis, cmap=cmap, figsize=figsize)\n        if y is not None: y.show(ax=ax, **kwargs)\n        if title is not None: ax.set_title(title)","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:07.11156Z","start_time":"2019-10-18T13:41:07.105636Z"},"trusted":true},"cell_type":"code","source":"def open_dicom_image(fn:PathOrStr, cls:type=DicomImage, \n                     after_open:Callable=None, expand=True,\n                     clamp_min=-1024, clamp_max=1024)->Image:\n    \"Return `Image` object created from image in file `fn`.\"\n    \n    ds = pydicom.dcmread(fn) # open dicom image as dicom dataset\n    img = ds.pixel_array # get pixel data as np array\n    \n    # Convert to Hounsfield units (HU) and rescale and set intercept.\n    # In this setup we only take a look at the values between -1024 and 1024.\n    # Values below will be set to -1024, values above to 1024\n    resc_img = img * ds.RescaleSlope + ds.RescaleIntercept\n    resc_img[resc_img < -1024] = clamp_min # Clamp to minimum value   \n    resc_img[resc_img > 1024] = clamp_max # Clamp to maximum value\n    resc_img = (resc_img - clamp_min) / (clamp_max - clamp_min) # rescale to range from 0 to 1\n    \n    if after_open: resc_img = after_open(resc_img)\n    \n    resc_img = torch.from_numpy(resc_img.astype(np.float32, copy=False))\n    \n    x = resc_img.view(1,*resc_img.shape)\n    \n    if expand: x = x.expand(3,*x.shape[-2:])\n    # x.shape[-2:] is needed because not everything is of size 512x512!\n    # expand is memory efficient: https://stackoverflow.com/questions/44593141/stacking-copies-of-an-array-a-torch-tensor-efficiently\n                \n    return cls(x)","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:07.098949Z","start_time":"2019-10-18T13:41:07.094575Z"},"trusted":true},"cell_type":"code","source":"def show_dicom_image(img:Image, ax:plt.Axes=None, figsize:tuple=(3,3), hide_axis:bool=True, cmap:str=cmap,\n                alpha:float=None, **kwargs)->plt.Axes:\n    \"Display `DicomImage` in the notebook.\"\n    if ax is None: fig,ax = plt.subplots(figsize=figsize)\n    ax.imshow(image2np(img.data.data[0].unsqueeze(0)), # We show only one channel!\n              cmap=cmap, alpha=alpha, **kwargs)\n    if hide_axis: ax.axis('off')\n    return ax","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:07.118912Z","start_time":"2019-10-18T13:41:07.112494Z"},"trusted":true},"cell_type":"code","source":"class DicomImageList(ImageList):\n    def __init__(self, *args, after_open:Callable=None, **kwargs):\n        super().__init__(*args, **kwargs)\n        self.after_open = after_open\n        self.c,self.sizes = 1,{}\n        \n    def open(self, fn):\n        \"Open image in `fn`, subclass and overwrite for custom behavior.\"\n        return open_dicom_image(fn, after_open=self.after_open)\n    \n    # based on https://github.com/radekosmulski/rsna-intracranial/blob/master/03_train_basic_model.ipynb\n    def label_from_df(self, cols:IntsOrStrs=1, label_cls:Callable=None, **kwargs):\n        \"Label `self.items` from the values in `cols` in `self.inner_df`.\"\n        self.inner_df.labels.fillna('', inplace=True)\n        labels = self.inner_df.iloc[:,df_names_to_idx(cols, self.inner_df)]\n        assert labels.isna().sum().sum() == 0, f\"You have NaN values in column(s) {cols} of your dataframe, please fix it.\"\n        if is_listy(cols) and len(cols) > 1 and (label_cls is None or label_cls == MultiCategoryList):\n            new_kwargs,label_cls = dict(one_hot=True, classes= cols),MultiCategoryList\n            kwargs = {**new_kwargs, **kwargs}\n        return self._label_from_list(_maybe_squeeze(labels), label_cls=label_cls, **kwargs)\n    \n    def show_xys(self, xs, ys, imgsize:int=4, figsize:Optional[Tuple[int,int]]=None, **kwargs):\n        \"Show the `xs` (inputs) and `ys` (targets) on a figure of `figsize`.\"\n        rows = int(np.ceil(math.sqrt(len(xs))))\n        axs = subplots(rows, rows, imgsize=imgsize, figsize=figsize)\n        for x,y,ax in zip(xs, ys, axs.flatten()): x.show(ax=ax, y=y, **kwargs)\n        for ax in axs.flatten()[len(xs):]: ax.axis('off')\n        plt.tight_layout()\n        \n    def reconstruct(self, t): return DicomImage(t)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Test open_dicom_image"},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:07.125489Z","start_time":"2019-10-18T13:41:07.119816Z"},"trusted":true},"cell_type":"code","source":"fn = '../input/rsna-intracranial-hemorrhage-detection/stage_1_train_images/ID_52c9913b1.dcm'","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:07.139509Z","start_time":"2019-10-18T13:41:07.126397Z"},"trusted":true},"cell_type":"code","source":"img = open_dicom_image(fn)","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:07.14245Z","start_time":"2019-10-18T13:41:07.140337Z"},"trusted":true},"cell_type":"code","source":"img.data.shape","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:07.208422Z","start_time":"2019-10-18T13:41:07.143209Z"},"trusted":true},"cell_type":"code","source":"img","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:07.222519Z","start_time":"2019-10-18T13:41:07.209253Z"},"trusted":true},"cell_type":"code","source":"# see if resize works properly\nimg.resize(128)","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:07.226224Z","start_time":"2019-10-18T13:41:07.223351Z"},"trusted":true},"cell_type":"code","source":"# check min & max\nimg.data.min(), img.data.max()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Databunch"},{"metadata":{},"cell_type":"markdown","source":"Now we can create our Databunch based on the `DicomImageList`."},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:07.229395Z","start_time":"2019-10-18T13:41:07.227041Z"},"trusted":true},"cell_type":"code","source":"# set image and batch size\nsz, bs = 256, 32","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:09.385402Z","start_time":"2019-10-18T13:41:07.230167Z"},"trusted":true},"cell_type":"code","source":"data = (DicomImageList.from_df(df_train,\n                               path='../input/rsna-intracranial-hemorrhage-detection',\n                               folder='stage_1_train_images')\n        .split_none()\n        .label_from_df(cols=-1, label_delim=' ')\n        .transform(size=sz)\n        .databunch(bs=bs))","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:09.388085Z","start_time":"2019-10-18T13:41:09.386318Z"},"trusted":true},"cell_type":"code","source":"# these dummy statistics are only calculated on three random batches\nstats_dicom = [torch.tensor(0.2192), torch.tensor(0.2775)]","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:09.417881Z","start_time":"2019-10-18T13:41:09.38957Z"},"trusted":true},"cell_type":"code","source":"data.normalize(stats_dicom)","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:09.420548Z","start_time":"2019-10-18T13:41:09.418842Z"},"trusted":true},"cell_type":"code","source":"# verify that we have our 6 classes\nassert data.c == 6","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now we are ready to have a look at the data!"},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:13.366755Z","start_time":"2019-10-18T13:41:09.42139Z"},"trusted":true},"cell_type":"code","source":"data.show_batch()","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:15.332754Z","start_time":"2019-10-18T13:41:13.367793Z"},"trusted":true},"cell_type":"code","source":"data.show_batch()","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-10-18T13:41:17.496146Z","start_time":"2019-10-18T13:41:15.333781Z"},"trusted":true},"cell_type":"code","source":"data.show_batch()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Feel free to improve and extend this kernel! Please comment if you have suggestion, ideas, have found errors, etc. I am happy if the kernel is useful to the community! :-D"}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.7.4"},"toc":{"base_numbering":1,"nav_menu":{},"number_sections":true,"sideBar":true,"skip_h1_title":false,"title_cell":"Table of Contents","title_sidebar":"Contents","toc_cell":false,"toc_position":{"height":"calc(100% - 180px)","left":"10px","top":"150px","width":"215px"},"toc_section_display":true,"toc_window_display":true}},"nbformat":4,"nbformat_minor":1}