{"cells":[{"metadata":{},"cell_type":"markdown","source":"This is the jhoward notebook saving to zip file to try to overcome the Kaggle output limit but it does not complete with the 5Gb data limit.\nI had to tweak a bit to have it working on Kaggle kernel :\n* save images as TIFF and not PNG to avoid exception of mode I;16 in PNG writer\n* save to zipfile through io.BytesIO\n"},{"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-png).)\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 16 bit PNG 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":{"trusted":true},"cell_type":"code","source":"from zipfile import ZipFile ","execution_count":null,"outputs":[]},{"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'\ntorch.set_num_threads(1)\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 16 bit images"},{"metadata":{},"cell_type":"markdown","source":"We will save our smaller cropped files as 16 bit images (formerly PNG, here TIFF). It's a really convenient format since you can open it with standard libraries. And 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 the 16 bit 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        #print('read', fn)\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)\npath_dest = path_dest/'train_png'\n# NB: Use bs=512 for GPUs with <16GB RAM\n#bs=1024\nbs=512\n\ndsrc = DataSource(fns, [[dcm_tfm],[os.path.basename]])\n# it worked with only 1 worker for the moment\ndl = TfmdDL(dsrc, bs=bs, num_workers=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\ndef save_img(x, buf, bins, compress_level, format=\"TIFF\"):\n    im = Image.fromarray(x.to_uint16(bins))\n    if compress_level: im.save(buf, compress_level=compress_level, format=format)\n    else: im.save(buf, compress_type=Image.RLE, format=format)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We'll need a way to save a file as 16-bit PNG - here it is! Note that we need to pass in `bins`, since it will use frequency-histogram normalization automatically with this bins."},{"metadata":{"trusted":true},"cell_type":"code","source":"def dest_fname_dest(fname, dest, suffix): return str(dest/Path(fname).with_suffix(suffix))\n    \ndef save_cropped_img(o, dest):\n    def dest_fname(fname): return dest_fname_dest(fname, dest, suffix='.tiff')\n    fname,px,zip_file = o\n    buf = io.BytesIO()\n    save_img(px, buf, bins=bins, compress_level=9, format='TIFF')\n    zip_file.writestr(dest_fname_dest(fname, dest), buf.getvalue())\n","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, dest, n_workers=4):\n    zip_files = []\n    for worker in range(n_workers):\n        zip_files.append(ZipFile(dest/'result_dataset_{}.zip'.format(worker), 'w'))\n    pxs = pxs.cuda()\n    masks = pxs.mask_from_blur(dicom_windows.brain)\n    bbs = mask2bbox(masks)\n    gs = crop_resize(pxs, bbs, 256).cpu().squeeze()\n    #g = partial(save_cropped_img, dest=dest)\n    parallel(save_cropped_img, zip(fnames, gs, zip_files), n_workers=n_workers, progress=False, dest=dest)\n    \ndef process_batch_no_parallel(pxs, fnames, dest, zip_file):\n    pxs = pxs.cuda()\n    masks = pxs.mask_from_blur(dicom_windows.brain)\n    bbs = mask2bbox(masks)\n    gs = crop_resize(pxs, bbs, 256).cpu().squeeze()\n    for fname, px in zip(fnames, gs):\n        buf = io.BytesIO()\n        save_img(px, buf, bins=bins, compress_level=9, format='TIFF')\n        zip_file.writestr(dest_fname_dest(fname, dest, suffix='.tiff'), buf.getvalue())\n","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-png) 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":"#!rm -f core\ndest = path_dest\ndest.mkdir(exist_ok=True)\n# , force_zip64=True\nwith ZipFile(dest/'result_dataset.zip', 'w' ) as zip_file:\n    for b in progress_bar(dl): process_batch_no_parallel(*b, dest=dest, zip_file=zip_file)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":" ## List of produced files and size"},{"metadata":{"trusted":true},"cell_type":"code","source":"dest = path_dest\ndest.mkdir(exist_ok=True)\n#ZipFile(dest/'result_dataset.zip', 'w')\nfor path in  ['.', str(dest)]:\n    for fname in os.listdir(path):\n        full_fname = os.path.join(path, fname)\n        print(full_fname, os.stat(full_fname).st_size)","execution_count":null,"outputs":[]},{"metadata":{"scrolled":false,"trusted":true},"cell_type":"code","source":"# Uncomment this to test and time a single batch\n\n# %time process_batch(*dl.one_batch(), n_workers=4)","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#     pngfn = dest/Path(fn).with_suffix('.png').name\n#     a = pngfn.png16read()\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}