{"cells":[{"metadata":{},"cell_type":"markdown","source":"Let's try to clean up the dataset as best as we can, and create something we can do some model prototyping with! For prototyping we'll want something that's small, and ideally fits within Kaggle's 20GB limit so we can upload it as a dataset for others to easily work with. (If you're just looking for the dataset created with this notebook, you can [find it here](https://www.kaggle.com/jhoward/rsna-hemorrhage-jpg).)\n\nAfter cleaning up the data and create a prototyping dataset, we will use it to train and submit a model. You can see the followup notebook that uses this data here: [From Prototyping To Submission](https://www.kaggle.com/jhoward/from-prototyping-to-submission-fastai)\n\nHere are the issues that we will address.\n\n1. Fix images with incorrect RescaleIntercept\n1. Remove some images if they have little useful information (e.g. they don't actually contain brain tissue)\n1. Resample this dataset to 2/1 split of with/without haemorrhage, so we have a smaller dataset for quick prototyping\n1. Crop the images to just contain the brain, and save the size of the crop in case it's important\n1. Do histogram rescaling and then save JPEG 256x256 px images\n\nWe'll be using the fastai.medical.imaging library here - for more information about this see the notebook [Some DICOM gotchas to be aware of](https://www.kaggle.com/jhoward/some-dicom-gotchas-to-be-aware-of-fastai). We'll also use the same basic setup that's in the notebook."},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"cell_type":"code","source":"!pip install torch torchvision feather-format kornia pyarrow --upgrade   > /dev/null\n!pip install git+https://github.com/fastai/fastai_dev                    > /dev/null","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"cell_type":"code","source":"from fastai2.basics           import *\nfrom fastai2.vision.all       import *\nfrom fastai2.medical.imaging  import *\n\nnp.set_printoptions(linewidth=120)\nmatplotlib.rcParams['image.cmap'] = 'bone'\nset_seed(42)\nset_num_threads(1)\n\npath = Path('../input/rsna-intracranial-hemorrhage-detection/')\npath_trn = path/'stage_1_train_images'\npath_tst = path/'stage_1_test_images'\npath_dest = Path()\npath_dest.mkdir(exist_ok=True)\n\npath_inp = Path('../input')","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"cell_type":"code","source":"path_df = path_inp/'creating-a-metadata-dataframe'\ndf_lbls = pd.read_feather(path_df/'labels.fth')\ndf_tst = pd.read_feather(path_df/'df_tst.fth')\ndf_trn = pd.read_feather(path_df/'df_trn.fth').dropna(subset=['img_pct_window'])\ncomb = df_trn.join(df_lbls.set_index('ID'), 'SOPInstanceUID')","execution_count":null,"outputs":[]},{"metadata":{"heading_collapsed":true},"cell_type":"markdown","source":"# Fix incorrect RescaleIntercept"},{"metadata":{"hidden":true},"cell_type":"markdown","source":"In an earlier notebook (\"[Some DICOM gotchas to be aware of](https://www.kaggle.com/jhoward/some-dicom-gotchas-to-be-aware-of-fastai)\") we saw that there are some images which seem to have incorrect values for RescaleIntercept. In addition, in [this comment](https://www.kaggle.com/jhoward/some-dicom-gotchas-to-be-aware-of-fastai/comments#646349) Malcolm McLean pointed out that the values in these images seem wrong. Let's see if we can figure out what's happening, and fix them! Here are the data subsets we created in the previous notebook:"},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","hidden":true,"trusted":true},"cell_type":"code","source":"repr_flds = ['BitsStored','PixelRepresentation']\ndf1 = comb.query('(BitsStored==12) & (PixelRepresentation==0)')\ndf2 = comb.query('(BitsStored==12) & (PixelRepresentation==1)')\ndf3 = comb.query('BitsStored==16')\ndfs = L(df1,df2,df3)","execution_count":null,"outputs":[]},{"metadata":{"hidden":true},"cell_type":"markdown","source":"The problematic images are those in `df1`, which don't have the expected `RescaleIntercept` of `-1024` or similar. We'll grab that subset, and have a look at a few of them"},{"metadata":{"hidden":true,"trusted":true},"cell_type":"code","source":"def df2dcm(df): return L(Path(o).dcmread() for o in df.fname.values)","execution_count":null,"outputs":[]},{"metadata":{"hidden":true,"scrolled":true,"trusted":true},"cell_type":"code","source":"df_iffy = df1[df1.RescaleIntercept>-100]\ndcms = df2dcm(df_iffy)\n\n_,axs = subplots(2,4, imsize=3)\nfor i,ax in enumerate(axs.flat): dcms[i].show(ax=ax)","execution_count":null,"outputs":[]},{"metadata":{"hidden":true},"cell_type":"markdown","source":"Oh dear, they don't look good at all! Let's pick one for further analysis, and start with looking at the histogram of pixel values:"},{"metadata":{"hidden":true,"trusted":true},"cell_type":"code","source":"dcm = dcms[2]\nd = dcm.pixel_array\nplt.hist(d.flatten());","execution_count":null,"outputs":[]},{"metadata":{"hidden":true},"cell_type":"markdown","source":"Normally the mode for unsigned data images is zero, since they are the background pixels, as you see here:"},{"metadata":{"hidden":true,"trusted":true},"cell_type":"code","source":"d1 = df2dcm(df1.iloc[[0]])[0].pixel_array\nplt.hist(d1.flatten());","execution_count":null,"outputs":[]},{"metadata":{"hidden":true},"cell_type":"markdown","source":"Instead, our mode is:"},{"metadata":{"hidden":true,"trusted":true},"cell_type":"code","source":"scipy.stats.mode(d.flatten()).mode[0]","execution_count":null,"outputs":[]},{"metadata":{"hidden":true},"cell_type":"markdown","source":"My guess is that what happened in the \"iffy\" images is that they were actually signed data, but were treated as unsigned. If that's the case, the a value of `-1000` or `-1024` (the usual values for background pixels in signed data images) will have wrapped around to `4096-1000=3096`. So we'll need to shift everything up by `1000`, then move the values larger than `2048` back to where they should have been."},{"metadata":{"hidden":true,"trusted":true},"cell_type":"code","source":"d += 1000\n\npx_mode = scipy.stats.mode(d.flatten()).mode[0]\nd[d>=px_mode] = d[d>=px_mode] - px_mode\ndcm.PixelData = d.tobytes()\ndcm.RescaleIntercept = -1000","execution_count":null,"outputs":[]},{"metadata":{"hidden":true,"trusted":true},"cell_type":"code","source":"plt.hist(dcm.pixel_array.flatten());","execution_count":null,"outputs":[]},{"metadata":{"hidden":true},"cell_type":"markdown","source":"Let's see if that helped."},{"metadata":{"hidden":true,"trusted":true},"cell_type":"code","source":"_,axs = subplots(1,2)\ndcm.show(ax=axs[0]);   dcm.show(dicom_windows.brain, ax=axs[1])","execution_count":null,"outputs":[]},{"metadata":{"hidden":true},"cell_type":"markdown","source":"That looks pretty much perfect! We'll put that into a function that we can use to fix all our problematic images."},{"metadata":{"hidden":true,"trusted":true},"cell_type":"code","source":"def fix_pxrepr(dcm):\n    if dcm.PixelRepresentation != 0 or dcm.RescaleIntercept<-100: return\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","execution_count":null,"outputs":[]},{"metadata":{"hidden":true},"cell_type":"markdown","source":"Let's see if they all clean up so nicely."},{"metadata":{"hidden":true,"trusted":true},"cell_type":"code","source":"dcms = df2dcm(df_iffy)\ndcms.map(fix_pxrepr)\n\n_,axs = subplots(2,4, imsize=3)\nfor i,ax in enumerate(axs.flat): dcms[i].show(ax=ax)","execution_count":null,"outputs":[]},{"metadata":{"heading_collapsed":true},"cell_type":"markdown","source":"# Remove useless images"},{"metadata":{"hidden":true},"cell_type":"markdown","source":"Our goal here is to create a small, fast, convenient dataset for rapid prototyping. So let's get rid of images that don't provide much useful information, such as those with very little actual brain tissue in them. Brain tissue is in the region `(0,80)`. Let's find out how many pixels in this region are in each image. When we created the metadata data frame, we got a `img_pct_window` column included which has the % of pixels in the brain window."},{"metadata":{"hidden":true,"trusted":true},"cell_type":"code","source":"df_iffy.img_pct_window[:10].values","execution_count":null,"outputs":[]},{"metadata":{"hidden":true},"cell_type":"markdown","source":"We see that the first image contains nearly no brain tissue. It seems unlikely that images like this will have noticable haemorrhages. Let's test this hypothesis."},{"metadata":{"hidden":true,"trusted":true},"cell_type":"code","source":"plt.hist(comb.img_pct_window,40);","execution_count":null,"outputs":[]},{"metadata":{"hidden":true},"cell_type":"markdown","source":"There are a *lot* of images with nearly no brain tissue in them - presumably they're the slices above and below the brain. Let's see if they have any labels:"},{"metadata":{"hidden":true,"scrolled":false,"trusted":true},"cell_type":"code","source":"comb = comb.assign(pct_cut = pd.cut(comb.img_pct_window, [0,0.02,0.05,0.1,0.2,0.3,1]))\ncomb.pivot_table(values='any', index='pct_cut', aggfunc=['sum','count']).T","execution_count":null,"outputs":[]},{"metadata":{"hidden":true},"cell_type":"markdown","source":"We can see that, as expected, the images with little brain tissue (<2% of pixels) have almost no labels. So let's remove them. (Interestingly, we can also see a strong relationship between these two variables.)"},{"metadata":{"hidden":true,"trusted":true},"cell_type":"code","source":"comb.drop(comb.query('img_pct_window<0.02').index, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"heading_collapsed":true},"cell_type":"markdown","source":"# Resample to 2/3 split"},{"metadata":{"hidden":true},"cell_type":"markdown","source":"We will keep every row with a label:"},{"metadata":{"hidden":true,"trusted":true},"cell_type":"code","source":"df_lbl = comb.query('any==True')\nn_lbl = len(df_lbl)\nn_lbl","execution_count":null,"outputs":[]},{"metadata":{"hidden":true},"cell_type":"markdown","source":"...and we'll keep half that number of images without a label, which should keep the resultant size under Kaggle's 20GB dataset limit:"},{"metadata":{"hidden":true,"trusted":true},"cell_type":"code","source":"df_nonlbl = comb.query('any==False').sample(n_lbl//2)\nlen(df_nonlbl)","execution_count":null,"outputs":[]},{"metadata":{"hidden":true},"cell_type":"markdown","source":"Let's put them altogether and see how many we have."},{"metadata":{"hidden":true,"trusted":true},"cell_type":"code","source":"comb = pd.concat([df_lbl,df_nonlbl])\nlen(comb)","execution_count":null,"outputs":[]},{"metadata":{"heading_collapsed":true},"cell_type":"markdown","source":"# Crop to just brain area"},{"metadata":{"hidden":true},"cell_type":"markdown","source":"To create a smaller and faster dataset, we'll need smaller images. So let's make sure they contain the important information, by cropping out the non-brain area. To do so, we start with an image like this:"},{"metadata":{"hidden":true,"trusted":true},"cell_type":"code","source":"dcm = Path(dcms[3].filename).dcmread()\nfix_pxrepr(dcm)","execution_count":null,"outputs":[]},{"metadata":{"hidden":true,"trusted":true},"cell_type":"code","source":"px = dcm.windowed(*dicom_windows.brain)\nshow_image(px);","execution_count":null,"outputs":[]},{"metadata":{"hidden":true},"cell_type":"markdown","source":"...then blur it, to remove the small and thin areas:"},{"metadata":{"hidden":true,"trusted":true},"cell_type":"code","source":"blurred = gauss_blur2d(px, 100)\nshow_image(blurred);","execution_count":null,"outputs":[]},{"metadata":{"hidden":true},"cell_type":"markdown","source":"...and just select the areas that are bright in this picture:"},{"metadata":{"hidden":true,"trusted":true},"cell_type":"code","source":"show_image(blurred>0.3);","execution_count":null,"outputs":[]},{"metadata":{"hidden":true},"cell_type":"markdown","source":"We can use `fastai`'s `mask_from_blur` method to do this for us. We'll overlay the results on a few images to see if it looks OK:"},{"metadata":{"hidden":true,"trusted":true},"cell_type":"code","source":"_,axs = subplots(1,4, imsize=3)\nfor i,ax in enumerate(axs.flat):\n    dcms[i].show(dicom_windows.brain, ax=ax)\n    show_image(dcms[i].mask_from_blur(dicom_windows.brain), cmap=plt.cm.Reds, alpha=0.6, ax=ax)","execution_count":null,"outputs":[]},{"metadata":{"hidden":true},"cell_type":"markdown","source":"It's not perfect, but it'll do for our prototyping purposes. Now we need something that finds the extreme pixels. That turns out to be fairly simple:"},{"metadata":{"hidden":true,"trusted":true},"cell_type":"code","source":"def pad_square(x):\n    r,c = x.shape\n    d = (c-r)/2\n    pl,pr,pt,pb = 0,0,0,0\n    if d>0: pt,pd = int(math.floor( d)),int(math.ceil( d))        \n    else:   pl,pr = int(math.floor(-d)),int(math.ceil(-d))\n    return np.pad(x, ((pt,pb),(pl,pr)), 'minimum')\n\ndef crop_mask(x):\n    mask = x.mask_from_blur(dicom_windows.brain)\n    bb = mask2bbox(mask)\n    if bb is None: return\n    lo,hi = bb\n    cropped = x.pixel_array[lo[0]:hi[0],lo[1]:hi[1]]\n    x.pixel_array = pad_square(cropped)","execution_count":null,"outputs":[]},{"metadata":{"hidden":true,"trusted":true},"cell_type":"code","source":"_,axs = subplots(1,2)\ndcm.show(ax=axs[0])\ncrop_mask(dcm)\ndcm.show(ax=axs[1]);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Save JPEG images"},{"metadata":{},"cell_type":"markdown","source":"We will save our smaller cropped files as JPEGs. All our metadata is already in data frames, so there's no need to store it inside the image files too, like DICOM does.\n\nWhen we save the files, we'll first use frequency-binned histogram normalization, so that one of the channels in the files will include a nice distribution of pixels, without needing any further processing. See this notebook to learn how this works: \"[DON'T see like a radiologist!](https://www.kaggle.com/jhoward/don-t-see-like-a-radiologist-fastai)\"."},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"cell_type":"code","source":"htypes = 'any','epidural','intraparenchymal','intraventricular','subarachnoid','subdural'\n\ndef get_samples(df):\n    recs = [df.query(f'{c}==1').sample() for c in htypes]\n    recs.append(df.query('any==0').sample())\n    return pd.concat(recs).fname.values\n\nsample_fns = concat(*dfs.map(get_samples))\nsample_dcms = tuple(Path(o).dcmread().scaled_px for o in sample_fns)\nsamples = torch.stack(sample_dcms)\nbins = samples.freqhist_bins()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We'll also save those bins, since we'll need them for processing the full dataset when we use it later, and for the test set when it's time to submit."},{"metadata":{"trusted":true},"cell_type":"code","source":"(path_dest/'bins.pkl').save(bins)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here's the steps to read a fix a single file, ensuring it's the standard 512x512 size (nearly all are that size already, but we need them to be consistent in later processing). Also, if there are any broken files, we'll skip them, by raising fastai's `SkipItemException` (which means \"don't use this file in the `DataLoader`)."},{"metadata":{"trusted":true},"cell_type":"code","source":"def dcm_tfm(fn): \n    fn = Path(fn)\n    try:\n        x = fn.dcmread()\n        fix_pxrepr(x)\n    except Exception as e:\n        print(fn,e)\n        raise SkipItemException\n    if x.Rows != 512 or x.Columns != 512: x.zoom_to((512,512))\n    return x.scaled_px","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"If we're not careful, processing 300GB of input data could take a *really* long time! To make it super fast, we'll do it in 3 steps:\n\n- Create a multiprocess `DataLoader` that reads, fixes, and rescales the DICOMs, returning batches of size `bs`\n- Loop through each batch, moving it to the GPU, and then using fastai's GPU-optimized masking, cropping, and resizing functions on the whole batch at once\n- For each batch, save each image in it in parallel, using fastai's `parallel` function.\n\nThis is the first time that I've shown how to combine parallel and GPU processing for preprocessing and saving medical images (and might be the first time it's been shown anywhere, at least using something this flexible and concise!)\n\nLet's start by setting our params, creating our [DataSource](http://dev.fast.ai/data.core.html#DataSource), and then wrapping that with a multiprocess transformed data loader."},{"metadata":{"trusted":true},"cell_type":"code","source":"fns = list(comb.fname.values)\ndest = path_dest/'train_jpg'\ndest.mkdir(exist_ok=True)\n# NB: Use bs=512 or 1024 when running on GPU\nbs=4\n\ndsrc = DataSource(fns, [[dcm_tfm],[os.path.basename]])\ndl = TfmdDL(dsrc, bs=bs, num_workers=2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We'll need a way to save a file as jpg - here it is! Note that we need to pass in `bins`, since it will use frequency-histogram normalization automatically with these bins as one of its channels."},{"metadata":{"trusted":true},"cell_type":"code","source":"def dest_fname(fname): return dest/Path(fname).with_suffix('.jpg')\n\ndef save_cropped_jpg(o, dest):\n    fname,px = o\n    px.save_jpg(dest_fname(fname), dicom_windows.brain, dicom_windows.subdural, bins=bins)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Finally, we can write our function to do the compute-intensive masking, cropping, and resizing on the GPU, and then spin of the parallel processing for saving."},{"metadata":{"trusted":true},"cell_type":"code","source":"def process_batch(pxs, fnames, n_workers=4):\n    pxs = to_device(pxs)\n    masks = pxs.mask_from_blur(dicom_windows.brain)\n    bbs = mask2bbox(masks)\n    gs = crop_resize(pxs, bbs, 256).cpu().squeeze()\n    parallel(save_cropped_jpg, zip(fnames, gs), n_workers=n_workers, progress=False, dest=dest)","execution_count":null,"outputs":[]},{"metadata":{"scrolled":false,"trusted":true},"cell_type":"code","source":"# test and time a single batch. It's ~100x faster on a GPU!\n%time process_batch(*dl.one_batch(), n_workers=3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fn = dest.ls()[0]\nim = Image.open(fn)\nfn","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"show_images(tensor(im).permute(2,0,1), titles=['brain','subdural','normalized'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's run it! I'll leave it commented out in this kernel because Kaggle won't let us save this many files for output - I've uploaded [a dataset](https://www.kaggle.com/jhoward/rsna-hemorrhage-jpg) with the processed files that you can use, or just download this notebook yourself and run it on your machine. It only takes 20 minutes on my GPU!"},{"metadata":{"trusted":true},"cell_type":"code","source":"# dest.mkdir(exist_ok=True)\n# for b in progress_bar(dl): process_batch(*b, n_workers=8)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Uncomment this to view some processed images\n\n# for i,(ax,fn) in enumerate(zip(subplots(2,4)[1].flat,fns)):\n#     jpgfn = dest/Path(fn).with_suffix('.jpg').name\n#     a = jpgfn.jpg16read()\n#     show_image(a,ax=ax)","execution_count":null,"outputs":[]}],"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"}},"nbformat":4,"nbformat_minor":1}