{"cells":[{"metadata":{},"cell_type":"markdown","source":"In this notebook we're going to explore the important issue of *windowing* from a rather different direction to what you've probably seen before... and will show you a new way to handle DICOM pixel rescaling that might just give your models a boost!\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)."},{"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\n\nfrom fastai2.basics           import *\nfrom fastai2.medical.imaging  import *\n\nnp.set_printoptions(linewidth=120)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"path = Path('../input/rsna-intracranial-hemorrhage-detection/')\npath_trn = path/'stage_1_train_images'\npath_tst = path/'stage_1_test_images'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As discussed in the excellent kernel [See like a Radiologist with Systematic Windowing](https://www.kaggle.com/dcstang/see-like-a-radiologist-with-systematic-windowing) (from which I derived the somewhat cheeky title of this current notebook!), radiologists use *windowing* to increase the contrast of images across particular bands of *hounsfield units*.\n\nBut *why* do they do this? It's because the human visual system can only see 100 levels of gradation of a single color (in this case, white/grey/black) - and even fewer with the limitations of a computer display - but there are up to `2**16 == 65536` levels in a 16 bit DICOM image. So a human can't possibly see all those levels at once in a single greyscale image. To see this, let's look at the same image the *See like a radiologist* notebook studies, without any windows:"},{"metadata":{"scrolled":false,"trusted":true},"cell_type":"code","source":"fname = path_trn/'ID_9d9cc6b01.dcm'\ndcm = fname.dcmread()\ndcm.show(scale=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"On the other hand, with brain windowing, we can clearly see the details of the brain:"},{"metadata":{"trusted":true},"cell_type":"code","source":"dcm.show(scale=dicom_windows.brain)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**But**, do *computers* need to use windowing? Every kernel I've seen so far, and every person in this competition I've spoken to, are all using windowing. But... computers don't have the problem that their perception systems can only see 100 levels at a time! So, they don't need windowing!\n\nPerhaps, you're thinking, we need to squish the image down into 256 levels so we can put it in an 8-bit image file? Or, like [some folks have demonstrated](https://www.kaggle.com/reppic/gradient-sigmoid-windowing) we could get 3x more levels by using a 3 channel image with different windows in each channel?\n\nNope, we don't even need to do that! Remember: a neural network takes as input *floating point* data, not integer data. Floats (as used in deep learning libraries) use 32 bits, and give us a high level of precision, particuarly for numbers close to zero. (If you're not familiar with how floats work \"under the hood\", check out the \"Floating point arithmetic\" section of [Computational Linear Algebra for Coders](https://github.com/fastai/numerical-linear-algebra/blob/master/nbs/1.%20Why%20are%20we%20here.ipynb) for Dr Rachel Thomas.)\n\nIf we use floats, this does mean that we can't read, write, and manipulate our images with the popular PIL python package (which only fully supports 8 bit data). And we can't save it in formats like JPEG, since that's 8 bit too. But we can use OpenCV or fastai.vision, both of which can operate directly on floating point data."},{"metadata":{},"cell_type":"markdown","source":"# Rescaling floating point data"},{"metadata":{},"cell_type":"markdown","source":"Note, however, that this doesn't mean that we can ignore scaling entirely. As shown in detail in [Deep Learning from the Foundations](https://course.fast.ai/part2.html), and particularly [lesson 10](https://course.fast.ai/videos/?lesson=10&t=3641), having well-scaled inputs is really important to getting good results from your neural net training. That means we want to see a good mix of values throughout the range of our data - e.g. something having approximately a normal or uniform distribution. Unfortunately, the pixels in our DICOM above don't show that at all:"},{"metadata":{"trusted":true},"cell_type":"code","source":"px = dcm.scaled_px.flatten()\nplt.hist(px, bins=40);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We see a highly bimodal distribution, with lots of background pixels at around `-1000`, and brain tissue pixels at a bit over `0`. But that's OK, we can normalize this using a simple trick - we can scale our pixel values using a non-linear mapping designed to give us an equal number of pixels in each range. To do that, we first need to split the range of pixel values into groups, such that each group has around the same number of pixels. `fastai` has a method to do this:"},{"metadata":{"trusted":true},"cell_type":"code","source":"bins = px.freqhist_bins(20)\nprint(bins)\nplt.hist(px, bins=bins);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Because there's some pixel values (like `1024`) which appear many times, we can't get a perfectly uniform result, but it's pretty close. So now we just need a way to scale our pixels evenly using these bins. We can use a simple function which connects these bins with line segments, like this:"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(bins, torch.linspace(0,1,len(bins)));","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"`fastai.medical.imaging` can apply that mapping for you, like so:"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(dcm.hist_scaled(), cmap=plt.cm.bone);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In fact, this is the default way of displaying a DICOM in fastai:"},{"metadata":{"trusted":true},"cell_type":"code","source":"dcm.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Creating a normalized dataset"},{"metadata":{},"cell_type":"markdown","source":"We can't use this rescaling directly, for one very important reason: the mapping, the way we've defined it, will very from image to image. That means that the pixel values that represent brain tissue, for instance, will be different from image to image. This will be difficult for our neural net to handle. It would be like using different normalization mean and standard deviation for each image, when dealing with normalization in regular image neural nets.\n\nThe fix is to create a mapping that is appropriate for a wide range of images. Let's create a mapping from a set of data that represents the three main subgroups we saw in the [Some DICOM gotchas](https://www.kaggle.com/jhoward/some-dicom-gotchas-to-be-aware-of-fastai) notebook (see the section \"*Looking at metadata - BitsStored and PixelRepresentation*\":\n\nTo find the appropriate images, we can use the files we created in the notebook [Creating a metadata DataFrame](https://www.kaggle.com/jhoward/creating-a-metadata-dataframe-fastai)."},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"cell_type":"code","source":"path_inp = Path('../input')\npath_df = path_inp/'creating-a-metadata-dataframe-fastai'\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')\n\ncomb = df_trn.join(df_lbls.set_index('ID'), 'SOPInstanceUID')\nrepr_flds = ['BitsStored','PixelRepresentation']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"These are the 3 subsets we identified in that notebook:"},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"cell_type":"code","source":"df1 = 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":{},"cell_type":"markdown","source":"For each of these subsets, we'll grab a random image with each label, and one with no labels."},{"metadata":{"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))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now we will read each of those image files:"},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_dcms = L(Path(o).dcmread() for o in sample_fns)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"...and put them in a tensor:"},{"metadata":{"trusted":true},"cell_type":"code","source":"samples = torch.stack(tuple(sample_dcms.attrgot('scaled_px')))\nsamples.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"...and get the bins from this set of images:"},{"metadata":{"trusted":true},"cell_type":"code","source":"bins = samples.freqhist_bins()\nplt.plot(bins, torch.linspace(0,1,len(bins)));","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"If we pass these bins to `show`, then the image will scale with these bins:"},{"metadata":{"trusted":true},"cell_type":"code","source":"dcm.show(bins)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"More importantly, we can pass them to `hist_scaled_px` to get the scaled float tensor, which we can pass to our model:"},{"metadata":{"trusted":true},"cell_type":"code","source":"dcm.hist_scaled(bins)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Note that since this is approximately uniformally distributed between zero and one, it's doesn't have the mean zero, standard deviation one, distribution that our model will want:"},{"metadata":{"trusted":true},"cell_type":"code","source":"scaled_samples = torch.stack(tuple(o.hist_scaled(bins) for o in sample_dcms))\nscaled_samples.mean(),scaled_samples.std()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So all you have to do is to include the `hist_scaled_px` transformation in your model input pipeline, and normalize it to mean zero, standard deviation one. We'll see an example in a future notebook."},{"metadata":{},"cell_type":"markdown","source":"# Afterword: a better windowing for humans?"},{"metadata":{},"cell_type":"markdown","source":"Interestingly, we can use a \"rainbow colormap\" to fully utilize our computer's ability to display color, our visual system's ability to discern seven million different colors, and our well-distributed normalized pixel data to see an enormous amount of detail in a single image. The image we looked at earlier has a subdural hemorrhage, which you can clearly see as bright purple coloration around the outside of the brain, particulary in the bottom left."},{"metadata":{"trusted":true},"cell_type":"code","source":"dcm.show(cmap=plt.cm.gist_ncar, figsize=(6,6))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Obviously this color mapping could do with plenty of cleanup (for instance the background should all be black, not a psychodelic pattern!), but the general idea is hopefully clear. I've asked a few radiologists what they think of this approach, and many think it would be a pretty interesting experiment to try out... If any radiologists reading this give it a go, be sure to let me know."}],"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":4,"nbformat_minor":1}